1 Introduction

The Higgs boson, an elementary particle, has been researched by the Large Hadron Collider (LHC) as one of the primary scientific goals. Combining the updated data of the ATLAS [1] and CMS [2] Collaborations, now its measured mass is \(m_{h^0}=125.09\pm 0.24 \mathrm {GeV}\) [3], which represent that the Higgs mechanism is compellent. In the Standard Model (SM), the lepton-flavor number is conserved. However, the neutrino oscillation experiments [4,5,6,7,8,9,10,11,12] have convinced that neutrinos possess tiny masses and mix with each other. So the individual lepton numbers \(L_i=L_e,~L_\mu ,~L_\tau \) are not exact symmetries at the electroweak scale. Furthermore, the presentation of the GIM mechanism makes the charged lepton flavor violating(CFLV) processes in the SM very tiny [13,14,15], such as \(Br_{SM}(l_j\rightarrow l_i\gamma )\sim 10^{-55}\) [16]. Therefore, if we observe the CLFV processes in future experiments, it is an obvious evidence of new physics beyond the SM.

Studying the CLFV processes is an effective way to explore new physics beyond the SM. MEG Collaboration gives out the current experiment upper bound of the CLFV process \(\mu \rightarrow e \gamma \), which is \(Br(\mu \rightarrow e \gamma )< 5.7\times 10^{-13}\) at \(90\%\) confidence level [17]. \(Br(\tau \rightarrow e \gamma )< 3.3\times 10^{-8}\) and \(\;Br(\tau \rightarrow \mu \gamma )< 4.4\times 10^{-8}\) are also shown in Ref. [18]. SINDRUM II collaboration has updated the sensitivity of the \(\mu -e\) conversion rate in Au nuclei \(CR(\mu \rightarrow e:\;_{79}^{197}Au)<7\times 10^{-13}\) [19]. The current experiment upper bounds for the \(\tau \) decays \(Br(\tau \rightarrow 3e)<2.7\times 10^{-8}\) and \(Br(\tau \rightarrow 3\mu )<2.1\times 10^{-8}\) have been shown by Particle Data Group [18]. Furthermore, a direct research for the 125.1 GeV Higgs boson decays including the CLFV, \(h^0\rightarrow l_i l_j\), has been given out by the CMS Collaboration [20, 21] and ATLAS Collaboration [22]. We show the corresponding experiment upper bounds for processes \(h^0\rightarrow l_i l_j\) in Table 1. Physicists do more research on the CLFV processes for \(l_j\rightarrow l_i \gamma \) decays, \(\mu -e\) conversion rates in nuclei, the \(\tau \) decays and \(h^0\rightarrow l_i l_j\) decays in models beyond the SM [23,24,25,26,27,28]. In our previous work, we have studied \(l_j\rightarrow l_i \gamma \), muon conversion to electron in nuclei, the \(\tau \) decays and \(h^0\rightarrow l_i l_j\) processes in the \(\mu \nu \)SSM [29,30,31]. We also have discussed \(l_j\rightarrow l_i \gamma \) processes, the \(\tau \) decays and muon conversion to electron in nuclei in the BLMSSM [32, 33]. In this work, we study the processes \(l_j\rightarrow l_i \gamma \), \(\mu -e\) conversion in Au nuclei, the \(\tau \) decays and \(h^0\rightarrow l_i l_j\) in the extended BLMSSM, which is named as the EBLMSSM [34].

Table 1 Present experiment limits for 125 GeV Higgs decays \(h^0\rightarrow l_i l_j\) with the CLFV

Extending the MSSM with the introduced local gauged B and L, one obtains the so called BLMSSM [35,36,37,38]. In the BLMSSM, the exotic lepton masses are obtained from the Yukawa couplings with the two Higgs doublets \(H_u\) and \(H_d\). The values of these exotic lepton masses are around 100 GeV, which can well anastomose the current experiment bounds. However, with the development of high energy physics experiments, we may obtain the heavier experiment lower bounds of the exotic lepton masses in the near future, which makes the BLMSSM model do not exist. Therefore, two exotic Higgs superfields, the \(SU(2)_L\) singlets \(\Phi _{NL}\) and \(\varphi _{NL}\), are considered to be added in the BLMSSM. With the introduced superfields \(\Phi _{NL}\) and \(\varphi _{NL}\), the exotic leptons can turn heavy and should be unstable. This new model is named as the extended BLMSSM (EBLMSSM) [34]. In order to make the exotic leptons unstable, we add superfields Y and \(Y'\) in the EBLMSSM.

In the BLMSSM, the dark matter (DM) candidates include the lightest mass eigenstate of \(X, X^\prime \) mixing and a four-component spinor \(\tilde{X}\) composed by the superpartners of \(X,X^\prime \). In the EBLMSSM, the DM candidates not only include above terms presented in the BLMSSM, but also contain new terms due to the new introduced superfields of \(Y,Y'\). So the lighter mass eigenstates of \(Y,Y^\prime \) mixing and spinor \(\tilde{Y}\) are DM candidates [34, 39]. In Sect. 4.2 of our previous work [34], we suppose the lightest mass eigenstate of \(Y, Y^\prime \) mixing as a DM candidate, and calculate the relic density \(\Omega _Dh^2\). In the reasonable parameter space, \(\Omega _Dh^2\) of \(Y_1\) can match the experiment results well.

The Higgs boson \(h^0\) is produced chiefly from the gluon fusion \((gg\rightarrow h^0)\) at the LHC. The leading order (LO) contributions originate from the one-loop diagrams. In the BLMSSM, we have studied the \( h^0\rightarrow gg\) process in our previous work [40], and the virtual top quark loops play the dominate roles. The EBLMSSM results for \( h^0\rightarrow gg\) are same as those in BLMSSM, which have been discussed in Ref. [34]. Being different from BLMSSM, the exotic leptons in EBLMSSM are more heavy and the exotic sleptons of the 4th and 5th generations mix together to form a \(4\times 4\) mass matrices. The LO contributions for \(h^0\rightarrow \gamma \gamma \) originate from the one-loop diagrams. In the EBLMSSM [34], we have studied the decay \(h^0\rightarrow \gamma \gamma \) in detail. The processes \(h^0\rightarrow VV, V=(Z,W)\) also have been researched in this new model. Considering the constraints from the parameter space of these researches, we study the processes \(l_j\rightarrow l_i\gamma \), muon conversion to electron in Au nuclei, \(\tau \) decays and \(h^0\rightarrow l_il_j\) in this work.

The outline of this paper is organized as follows. In Sect. 2, we present the ingredients of the EBLMSSM by introducing its superpotential, the general soft SUSY-breaking terms, new corrected mass matrices and couplings which are different from those in the BLMSSM. In Sect. 3, we analyze the corresponding amplitudes and the branching ratios of rare CLFV processes \(l_j\rightarrow l_i \gamma \), the \(\tau \) decays and \(h^0\rightarrow l_i l_j\) decays. We also discuss the muon conversion to electron rates in nuclei. The numerical analysis is discussed in Sect. 4, and the conclusions are summarized in Sect. 5. The tedious formulae are collected in “Appendix”.

Table 2 The new introduced superfields in the EBLMSSM beyond BLMSSM

2 Introduction of the EBLMSSM

In the EBLMSSM, the local gauge group is \(SU(3)_{C}\otimes SU(2)_{L}\otimes U(1)_{Y}\otimes U(1)_{B}\otimes U(1)_{L}\) [34, 35, 41, 42]. We introduce the exotic Higgs superfields \(\Phi _{NL}\) and \(\varphi _{NL}\) with nonzero VEVs \(\upsilon _{NL}\) and \(\bar{\upsilon }_{NL}\) [43] to make the exotic leptons heavy. Accordingly, the superfields Y and \(Y'\) are introduced to avoid the heavy exotic leptons stable. In Table 2, we show the new introduced superfields in the EBLMSSM [34].

The corresponding superpotential of the EBLMSSM is shown here

$$\begin{aligned}&\mathcal{W}_{{EBLMSSM}}=\mathcal{W}_{{MSSM}}+\mathcal{W}_{B}+\mathcal{W}_{L}+\mathcal{W}_{X}+\mathcal{W}_{Y}, \nonumber \\&\mathcal{W}_{L}=\lambda _{L}\hat{L}_{4}\hat{L}_{5}^c\hat{\varphi }_{NL} \nonumber \\&\quad +\lambda _{E}\hat{E}_{4}^c\hat{E}_{5} \hat{\Phi }_{NL} \nonumber \\&\quad +\lambda _{NL}\hat{N}_{4}^c\hat{N}_{5}\hat{\Phi }_{NL} \nonumber \\&\quad +\mu _{NL}\hat{\Phi }_{NL}\hat{\varphi }_{NL} \nonumber \\&\quad +Y_{{e_4}}\hat{L}_{4}\hat{H}_{d}\hat{E}_{4}^c+Y_{{\nu _4}}\hat{L}_{4}\hat{H}_{u}\hat{N}_{4}^c \nonumber \\&\quad +Y_{{e_5}}\hat{L}_{5}^c\hat{H}_{u}\hat{E}_{5}+Y_{{\nu _5}}\hat{L}_{5}^c\hat{H}_{d}\hat{N}_{5}\nonumber \\&\quad +Y_{\nu }\hat{L}\hat{H}_{u}\hat{N}^c+\lambda _{{N^c}}\hat{N}^c\hat{N}^c\hat{\varphi }_{L} +\mu _{L}\hat{\Phi }_{L}\hat{\varphi }_{L}\;, \nonumber \\&\mathcal{W}_{Y}=\lambda _4\hat{L}\hat{L}_{5}^c\hat{Y}+\lambda _5\hat{N}^c\hat{N}_{5}\hat{Y}^\prime +\lambda _6\hat{E}^c\hat{E}_{5}\hat{Y}^\prime +\mu _{Y}\hat{Y}\hat{Y}^\prime \;.\nonumber \\ \end{aligned}$$
(1)

\(\mathcal{W}_{{MSSM}}\) is the superpotential of the MSSM. \(\mathcal{W}_{B}\) and \(\mathcal{W}_{X}\) are same as the terms in the BLMSSM [40, 44]. The new terms \(\lambda _{L}\hat{L}_{4}\hat{L}_{5}^c\hat{\varphi }_{NL}+\lambda _{E}\hat{E}_{4}^c\hat{E}_{5} \hat{\Phi }_{NL}+\lambda _{NL}\hat{N}_{4}^c\hat{N}_{5}\hat{\Phi }_{NL} +\mu _{NL}\hat{\Phi }_{NL}\hat{\varphi }_{NL}\) are added to \(\mathcal{W}_{L}\) based on the original BLMSSM. Comparing with the \(\mathcal{W}_{X}\) in the BLMSSM, \(\mathcal{W}_{Y}\) is introduced in the EBLMSSM, which includes the lepton-exotic lepton-Y coupling and lepton-exotic slepton-\(\tilde{Y}\) coupling. These new couplings can produce one-loop diagrams influencing the CLFV decays. These new couplings can also produce one-loop diagrams contributing to the lepton electric dipole moment(EDM) and lepton magnetic dipole moment(MDM), which will be discussed in our next work. With the 4th and 5th generation exotic sleptons mixing together, the \(h^0(Z)\)-exotic slepton-exotic slepton coupling is deduced in the EBLMSSM. In the EBLMSSM, the couplings for lepton–slepton–lepton neutralino, \(h^0(Z)\)-slepton–slepton, \(h^0(Z)\)-sneutrino–sneutrino and \(h^0(Z)\)-exotic lepton-exotic lepton also have new contributions to CLFV processes. In the whole, the new couplings in the EBLMSSM enrich the lepton physics in a certain degree.

In the EBLMSSM, \(\mathcal{W}_{Y}\) are the new terms in the superpotential. In \(\mathcal{W}_{Y}\), \(\lambda _4(\lambda _6)\) is the coupling coefficient of Y-lepton–exotic lepton and \(\tilde{Y}\)-slepton–exotic slepton couplings. We consider \(\lambda _4^2(\lambda _6^2)\) is a \(3\times 3\) matrix and has non-zero elements relating with the CLFV. In our following numerical analysis, we assume that \((\lambda _4^2)^{IJ}=(\lambda _6^2)^{IJ}=(Lm^2)^{IJ}\), I(J) represents the Ith (Jth) generation charged lepton. When \(I=J\), there is no CLFV, which has no contributions to our researched decay processes. So, only the non-diagonal elements \((Lm^2)^{IJ}(I\ne J)\) influence the numerical results of the CLFV processes. Therefore, we should take into account the effects from \(\mathcal{W}_{Y}\) in this work.

Based on the new introduced superfields \(\Phi _{NL},\varphi _{NL}, Y\) and \(Y'\) in the EBLMSSM, the soft breaking terms are given out

$$\begin{aligned}&\mathcal{L}_{{soft}}^{EBLMSSM}=\mathcal{L}_{{soft}}^{BLMSSM} -m_{{\Phi _{NL}}}^2\Phi _{NL}^*\Phi _{NL} \nonumber \\&\quad -m_{{\varphi _{NL}}}^2\varphi _{NL}^*\varphi _{NL} +(A_{{LL}}\lambda _{L}\tilde{L}_{4}\tilde{L}_{5}^c\varphi _{NL} \nonumber \\&\quad +A_{{LE}}\lambda _{E}\tilde{e}_{4}^c\tilde{e}_{5}\Phi _{NL} +A_{{LN}}\lambda _{NL}\tilde{\nu }_{4}^c\tilde{\nu }_{5}\Phi _{NL} \nonumber \\&\quad +B_{NL}\mu _{NL}\Phi _{NL}\varphi _{NL}+h.c.) \nonumber \\&\quad +(A_4\lambda _4\tilde{L}\tilde{L}_{5}^cY+A_5\lambda _5\tilde{N}^c\tilde{\nu }_{5}Y^\prime \nonumber \\&\quad +A_6\lambda _6\tilde{e}^c\tilde{e}_{5}Y^\prime +B_{Y}\mu _{Y}YY^\prime +h.c.). \end{aligned}$$
(2)

\(\mathcal{L}_{{soft}}^{BLMSSM}\) is the soft breaking terms of the BLMSSM discussed in our previous work [40, 44]. Here, corresponding to the \(SU(2)_L\) singlets \(\Phi _{NL}\) and \(\varphi _{NL}\), we obtain the nonzero VEVs \(\upsilon _{NL} \) and \(\bar{\upsilon }_{NL}\) respectively. Generally, the values of these two parameters are at TeV scale. The exotic Higgs \(\Phi _{NL}\) and \(\varphi _{NL}\) can be written as

$$\begin{aligned}&\Phi _{NL}={1\over \sqrt{2}}\Big (\upsilon _{NL}+\Phi _{NL}^0+iP_{NL}^0\Big ), \nonumber \\&\quad \varphi _{NL}={1\over \sqrt{2}}\Big (\bar{\upsilon }_{NL}+\varphi _{NL}^0+i\bar{P}_{NL}^0\Big ), \end{aligned}$$
(3)

where \(\tan \beta _{NL}=\bar{\upsilon }_{NL}/\upsilon _{NL}\) and \(v_{Nlt}=\sqrt{v_{NL}^2+\bar{v}_{NL}^2}\).

Comparing with the BLMSSM, the introduced superfields \(\Phi _{NL}\) and \(\varphi _{NL}\) in the EBLMSSM can give corrections to the mass matrices of the slepton, sneutrino, exotic lepton, exotic neutrino, exotic slepton, exotic sneutrino and lepton neutralino. However, the mass matrices of squark, exotic quark, exotic squark used in this work are same as those in the BLMSSM [40, 45]. We deduce the adjusted mass matrices in the EBLMSSM as follows.

2.1 The mass matrices of slepton and sneutrino in the EBLMSSM

In our previous work, we can easily obtain the slepton and sneutrino mass squared matrices of the BLMSSM [32]. Using the replacement \(\overline{\upsilon }^2_L-\upsilon ^2_L\rightarrow V_L^2\) (here \(V_L^2=\overline{\upsilon }^2_L-\upsilon ^2_L+\frac{3}{2}(\overline{\upsilon }^2_{NL}-\upsilon ^2_{NL})\)) for the BLMSSM results, we acquire the mass squared matrices of slepton and sneutrino in the EBLMSSM.

2.2 The mass matrices of exotic lepton and exotic neutrino in the EBLMSSM

The EBLMSSM exotic leptons masses are heavier than those in the BLMSSM due to the introduction of large parameters \(\upsilon _{NL}\) and \(\bar{\upsilon }_{NL}\). One can obtain the mass matrix of exotic lepton in the Lagrangian:

$$\begin{aligned}&-\mathcal{L}_{L'}^{mass}=\left( \begin{array}{ll}\bar{e}_{{4R}}^\prime&\bar{e}_{{5R}}^\prime \end{array}\right) \nonumber \\&\quad \left( \begin{array}{ll}-{1\over \sqrt{2}}\lambda _{L}\overline{\upsilon }_{NL}&{}{1\over \sqrt{2}}Y_{{e_5}}\upsilon _{u}\\ -{1\over \sqrt{2}}Y_{{e_4}}\upsilon _{d}&{}{1\over \sqrt{2}}\lambda _{E}\upsilon _{NL} \end{array}\right) \left( \begin{array}{l}e_{{4L}}^\prime \\ e_{{5L}}^\prime \end{array}\right) +h.c. \end{aligned}$$
(4)

Similarly, the mass matrix of the exotic neutrinos in the EBLMSSM can be given through the Lagrangian:

$$\begin{aligned}&-\mathcal{L}_{N'}^{mass}=\left( \begin{array}{ll}\bar{\nu }_{{4R}}^\prime&\bar{\nu }_{{5R}}^\prime \end{array}\right) \left( \begin{array}{ll}{1\over \sqrt{2}}\lambda _{L}\overline{\upsilon }_{NL}&{}-{1\over \sqrt{2}}Y_{{\nu _5}}\upsilon _{d}\\ {1\over \sqrt{2}}Y_{{\nu _4}}\upsilon _{u}&{}{1\over \sqrt{2}}\lambda _{NL}\upsilon _{NL} \end{array}\right) \nonumber \\&\left( \begin{array}{l}\nu _{{4L}}^\prime \\ \nu _{{5L}}^\prime \end{array}\right) +h.c. \end{aligned}$$
(5)

2.3 The mass matrices of exotic slepton and exotic sneutrino in the EBLMSSM

In the EBLMSSM, the exotic slepton of 4th generation and 5th generation mix together, and its mass matrix is \(4\times 4\), which is different from that in the BLMSSM. Using the superpotential in Eq. (1) and the soft breaking terms in Eq. (2), the mass squared matrix for exotic slepton can be obtained through Lagrangian:

$$\begin{aligned} -\mathcal{L}_{{\tilde{E}}}^{mass}=\tilde{E}^{\dag }\cdot \mathcal {M}^2_{\tilde{E}}\cdot \tilde{E}. \end{aligned}$$
(6)

With the base \(\tilde{E}^{T}=(\tilde{e}_4,\tilde{e}_4^{c*},\tilde{e}_5,\tilde{e}_5^{c*})\), we show the concrete elements of exotic slepton mass matrix \(\mathcal {M}^2_{\tilde{E}}\) in the following form

$$\begin{aligned}&\mathcal {M}^2_{\tilde{E}}(\tilde{e}_5^{c*}\tilde{e}_5^{c})= \lambda _L^2\frac{\bar{\upsilon }_{NL}^2}{2}+\frac{\upsilon _u^2}{2}|Y_{e_5}|^2+M^2_{\tilde{L}_5} \nonumber \\&\quad -\frac{g_1^2-g_2^2}{8}(\upsilon _d^2-\upsilon _u^2)-g_L^2(3+L_4)V_L^2, \nonumber \\&\mathcal {M}^2_{\tilde{E}}(\tilde{e}_5^{*}\tilde{e}_5)=\lambda _E^2\frac{\upsilon _{NL}^2}{2} +\frac{\upsilon _u^2}{2}|Y_{e_5}|^2+M^2_{\tilde{e}_5} \nonumber \\&\quad +\frac{g_1^2}{4}(\upsilon _d^2-\upsilon _u^2) +g_L^2(3+L_4)V_L^2, \nonumber \\&\mathcal {M}^2_{\tilde{E}}(\tilde{e}_4^{*}\tilde{e}_4)=\lambda _L^2\frac{\bar{\upsilon }_{NL}^2}{2} +\frac{g_1^2-g_2^2}{8}(\upsilon _d^2-\upsilon _u^2) \nonumber \\&\quad +\frac{\upsilon _d^2}{2}|Y_{e_4}|^2+M^2_{\tilde{L}_4} +g_L^2L_4V_L^2, \nonumber \\&\mathcal {M}^2_{\tilde{E}}(\tilde{e}_4^{c*}\tilde{e}_4^{c}) = \lambda _E^2\frac{\upsilon _{NL}^2}{2}-\frac{g_1^2}{4}(\upsilon _d^2-\upsilon _u^2) \nonumber \\&\quad +\frac{\upsilon _d^2}{2}|Y_{e_4}|^2+M^2_{\tilde{e}_4} -g_L^2L_4V_L^2, \nonumber \\&\mathcal {M}^2_{\tilde{E}}(\tilde{e}_4^{*}\tilde{e}_5) =\upsilon _dY_{e_4}^*\lambda _E\frac{\upsilon _{NL}}{2} +\lambda _LY_{e_5}\frac{\bar{\upsilon }_{NL}v_u}{2}, \nonumber \\&\quad \mathcal {M}^2_{\tilde{E}}(\tilde{e}_5\tilde{e}_5^{c})=\mu ^*\frac{\upsilon _d}{\sqrt{2}}Y_{e_5} +A_{e_5}Y_{e_5}\frac{\upsilon _u}{\sqrt{2}}, \nonumber \\&\mathcal {M}^2_{\tilde{E}}(\tilde{e}_4^{c}\tilde{e}_5)=\mu _{NL}^*\lambda _E \frac{\bar{\upsilon }_{NL}}{\sqrt{2}} -A_{LE}\lambda _E\frac{\upsilon _{NL}}{\sqrt{2}},\nonumber \\&~~\mathcal {M}^2_{\tilde{E}}(\tilde{e}_4\tilde{e}_5^{c}) = -\mu _{NL}^*\frac{\upsilon _{NL}}{\sqrt{2}}\lambda _L+A_{LL}\lambda _L\frac{\bar{\upsilon }_{NL}}{\sqrt{2}}, \nonumber \\&\mathcal {M}^2_{\tilde{E}}(\tilde{e}_4\tilde{e}_4^{c}) =\mu ^*\frac{\upsilon _u}{\sqrt{2}}Y_{e_4}+A_{e_4}Y_{e_4}\frac{\upsilon _d}{\sqrt{2}}, \nonumber \\&\quad \mathcal {M}^2_{\tilde{E}}(\tilde{e}_5^{c}\tilde{e}_4^{c*}) =-Y_{e_5}\lambda _E\frac{\upsilon _u\upsilon _{NL}}{2}-\lambda _LY_{e_4}^*\frac{\bar{\upsilon }_{NL}v_d}{2}. \end{aligned}$$
(7)

The matrix \(Z_{\tilde{E}}\) is used to rotate exotic slepton mass matrix to mass eigenstates, which is \(Z^{\dag }_{\tilde{E}}\mathcal {M}^2_{\tilde{E}} Z_{\tilde{E}}=diag(m^2_{\tilde{E}^1},m^2_{\tilde{E}^2},m^2_{\tilde{E}^3},m^2_{\tilde{E}^4})\).

In the same way, the exotic sneutrino mass squared matrix is also obtained through the Lagrangian:

$$\begin{aligned} -\mathcal{L}_{{\tilde{N}}}^{mass}=\tilde{N}^{\dag }\cdot \mathcal {M}^2_{\tilde{N}}\cdot \tilde{N}, \end{aligned}$$
(8)

where the corresponding elements of the matrix \(\mathcal {M}^2_{\tilde{N}}\) are

$$\begin{aligned}&\mathcal {M}^2_{\tilde{N}}(\tilde{\nu }_5^{c*}\tilde{\nu }_5^{c}) =\lambda _L^2\frac{\bar{\upsilon }_{NL}^2}{2} -\frac{g_1^2+g_2^2}{8}(\upsilon _d^2-\upsilon _u^2) \nonumber \\&\quad +\frac{\upsilon _d^2}{2}|Y_{\nu _5}|^2+M^2_{\tilde{L}_5} -g_L^2(3+L_4)V_L^2, \nonumber \\&\mathcal {M}^2_{\tilde{N}}(\tilde{\nu }_4^{*}\tilde{\nu }_4) =\lambda _L^2\frac{\bar{\upsilon }_{NL}^2}{2} +\frac{g_1^2+g_2^2}{8}(\upsilon _d^2-\upsilon _u^2) \nonumber \\&\quad +\frac{\upsilon _u^2}{2}|Y_{\nu _4}|^2+M^2_{\tilde{L}_4} +g_L^2L_4V_L^2, \mathcal {M}^2_{\tilde{N}}(\tilde{\nu }_5^{*}\tilde{\nu }_5) \nonumber \\&\quad = \lambda _{NL}^2\frac{\upsilon _{NL}^2}{2}+g_L^2(3+L_4)V_L^2 +\frac{\upsilon _d^2}{2}|Y_{\nu _5}|^2+M^2_{\tilde{\nu }_5}, \nonumber \\&\mathcal {M}^2_{\tilde{N}}(\tilde{\nu }_4^{c*}\tilde{\nu }_4^{c}) = \lambda _{NL}^2\frac{\upsilon _{NL}^2}{2}-g_L^2L_4V_L^2 +\frac{\upsilon _u^2}{2}|Y_{\nu _4}|^2+M^2_{\tilde{\nu }_4}, \nonumber \\&\mathcal {M}^2_{\tilde{N}}(\tilde{\nu }_5^{c}\tilde{\nu }_4^{c*}) = \lambda _{NL}Y_{\nu _5}\frac{\upsilon _{NL}\upsilon _d}{2} -\lambda _LY_{\nu _4}^*\frac{\bar{\upsilon }_{NL}\upsilon _u}{2}, \nonumber \\&\mathcal {M}^2_{\tilde{N}}(\tilde{\nu }_5\tilde{\nu }_5^{c}) = \mu ^*\frac{\upsilon _u}{\sqrt{2}}Y_{\nu _5}+A_{\nu _5}Y_{\nu _5}\frac{\upsilon _d}{\sqrt{2}}, \nonumber \\&\mathcal {M}^2_{\tilde{N}}(\tilde{\nu }_4^{c}\tilde{\nu }_5) =\mu _{NL}^*\lambda _{NL}\frac{\bar{\upsilon }_{NL}}{\sqrt{2}} -A_{LN}\lambda _N\frac{\upsilon _{NL}}{\sqrt{2}}, \nonumber \\&\mathcal {M}^2_{\tilde{N}}(\tilde{\nu }_4\tilde{\nu }_5^{c}) =\mu _{NL}^*\frac{\upsilon _{NL}}{\sqrt{2}} \lambda _L-A_{LL}\lambda _L\frac{\bar{\upsilon }_{NL}}{\sqrt{2}}, \nonumber \\&\mathcal {M}^2_{\tilde{N}}(\tilde{\nu }_4^{*}\tilde{\nu }_5) =\lambda _LY_{\nu _5}\frac{\bar{\upsilon }_{NL}\upsilon _d}{2} -\frac{\upsilon _u\upsilon _{NL}}{2}\lambda _{NL}Y_{\nu _4}^*, \nonumber \\&quad \mathcal {M}^2_{\tilde{N}}(\tilde{\nu }_4\tilde{\nu }_4^{c}) =\mu ^*\frac{\upsilon _d}{\sqrt{2}}Y_{\nu _4}+A_{\nu _4}Y_{\nu _4}\frac{\upsilon _u}{\sqrt{2}}. \end{aligned}$$
(9)

In the base \((\tilde{\nu }_4,\tilde{\nu }_4^{c*},\tilde{\nu }_5,\tilde{\nu }_5^{c*})\), we can diagonalize the mass squared matrix \(\mathcal {M}^2_{\tilde{N}}\) by \(Z_{\tilde{N}}\).

2.4 The lepton neutralino mass matrix in the EBLMSSM

In the EBLMSSM, \(\lambda _L\), the superpartner of the new lepton type gauge boson \(Z^\mu _L\), mixes with \((\psi _{\Phi _L},\psi _{\varphi _L},\psi _{\Phi _{NL}},\psi _{\varphi _{NL}})\) (the SUSY superpartners of the superfields (\(\Phi _{L},\varphi _{L},\Phi _{NL},\varphi _{NL}\))). So the lepton neutralino mass matrix is obtained in the base \((i\lambda _L,\psi _{\Phi _L},\psi _{\varphi _L},\psi _{\Phi _{NL}},\psi _{\varphi _{NL}})\),

$$\begin{aligned} \mathcal {M}_L=\left( \begin{array}{ccccc} 2M_L &{}2\upsilon _Lg_L &{}-2\bar{\upsilon }_Lg_L&{}3\upsilon _{NL}g_L &{}-3\bar{\upsilon }_{NL}g_L\\ 2\upsilon _Lg_L &{} 0 &{}-\mu _L&{} 0 &{} 0\\ -2\bar{\upsilon }_Lg_L&{}-\mu _L &{}0&{} 0 &{} 0\\ 3\upsilon _{NL}g_L &{} 0 &{} 0 &{} 0 &{} -\mu _{NL}\\ -3\bar{\upsilon }_{NL}g_L&{} 0&{}0&{}-\mu _{NL}&{}0 \end{array}\right) . \end{aligned}$$
(10)

The mass matrix \(\mathcal {M}_L\) can be diagonalized by the rotation matrix \(Z_{NL}\). Then, we can have

$$\begin{aligned}&i\lambda _L=Z_{NL}^{1i}K_{L_i}^0 ,~~~\psi _{\Phi _L}=Z_{NL}^{2i}K_{L_i}^0 ,~~~\psi _{\varphi _L}=Z_{NL}^{3i}K_{L_i}^0, \nonumber \\&\psi _{\Phi _{NL}}=Z_{NL}^{4i}K_{L_i}^0 ,~~~~~\psi _{\varphi _{NL}}=Z_{NL}^{5i}K_{L_i}^0. \end{aligned}$$
(11)

Here, \(X^0_{L_i}=(K_{L_i}^0,\bar{K}_{L_i}^0)^T\) represent the mass eigenstates of the lepton neutralino.

2.5 The superfields Y in the EBLMSSM

The scalar superfields Y and \(Y'\) mix. Adopting the unitary transformation,

$$\begin{aligned} \left( \begin{array}{c} Y_{1} \\ Y_{2}\\ \end{array}\right) =Z_{Y}^{\dag }\left( \begin{array}{c} Y \\ Y'^*\\ \end{array}\right) , \end{aligned}$$
(12)

the mass squared matrix for the superfield Y is deduced. With \(S_{Y}=g_{L}^2(2+L_{4})V_L^2\), the concrete form for the Y mass squared matrix is shown here

$$\begin{aligned} \mathcal{M}_Y^2=\left( \begin{array}{cc} |\mu _{Y}|^2+S_{Y} &{}-\mu _{Y}B_{Y} \\ -\mu ^*_{Y}B^*_{Y} &{} |\mu _{Y}|^2-S_{Y}\\ \end{array}\right) . \end{aligned}$$
(13)

The matrix \(Z_Y\) is used to diagonalize the matrix to the mass eigenstates:

$$\begin{aligned} Z^{\dag }_{Y}\left( \begin{array}{cc} |\mu _{Y}|^2+S_{Y} &{}-\mu _{Y}B_{Y} \\ -\mu ^*_{Y}B^*_{Y} &{} |\mu _{Y}|^2-S_{Y}\\ \end{array}\right) Z_{Y}=\left( \begin{array}{cc} m_{{Y_1}}^2 &{}0 \\ 0 &{} m_{{Y_2}}^2\\ \end{array}\right) . \end{aligned}$$
(14)

We suppose \(m_{{Y_1}}^2<m_{{Y_2}}^2\). The superpartners of Y and \(Y'\) form a four-component Dirac spinor \(\tilde{Y}\), and the mass term for superfield \(\tilde{Y}\) in the Lagrangian is given out

$$\begin{aligned}&-\mathcal {L}^{mass}_{\tilde{Y}}=\mu _Y\bar{\tilde{Y}}\tilde{Y} ,~~~~~~~~~~~~~~~~\tilde{Y} =\left( \begin{array}{c} \psi _{Y'} \\ \bar{\psi }_{Y}\\ \end{array}\right) . \end{aligned}$$
(15)
Fig. 1
figure 1

The triangle type diagrams for decays \(l_j\rightarrow l_i \gamma \)

In the EBLMSSM, the 4th and 5th generation exotic sleptons mix together. So the exotic slepton couplings in this new model are different from those in the BLMSSM. We deduce the \(h^0\)-exotic slepton-exotic slepton (\(h^0-\tilde{E}-\tilde{E}\)) coupling as follows

$$\begin{aligned} \mathcal {L}_{h^0\tilde{E}\tilde{E}}= & {} \sum _{i,j=1}^4\tilde{E}^{i*}\tilde{E}^{j}h^0 \nonumber \\&\times \Big [\Big (e^2\upsilon \sin \beta \frac{1-4s_W^2}{4s_W^2c_W^2}(Z_{\tilde{E}}^{4i*}Z_{\tilde{E}}^{4j} -Z_{\tilde{E}}^{1i*}Z_{\tilde{E}}^{1j}) \nonumber \\&-\frac{\mu ^*}{\sqrt{2}}Y_{e_4}Z_{\tilde{E}}^{2i*}Z_{\tilde{E}}^{1j} -\upsilon \sin \beta |Y_{e_5}|^2\delta _{ij}-\frac{A_{E_5}}{\sqrt{2}}Z_{\tilde{E}}^{4i*}Z_{\tilde{E}}^{3j} \nonumber \\&+\frac{1}{2}\lambda _LY_{e_5}Z_{\tilde{E}}^{3j}Z_{\tilde{E}}^{3i*}\bar{\upsilon }_{NL} -\frac{1}{2}Y_{e_5}^*Z_{\tilde{E}}^{4j}\lambda _EZ_{\tilde{E}}^{2i*}\upsilon _{NL}\Big )\cos \alpha \nonumber \\&- \Big (e^2\upsilon \cos \beta \frac{1-4s_W^2}{4s_W^2c_W^2}(Z_{\tilde{E}}^{1i*}Z_{\tilde{E}}^{1j}-Z_{\tilde{E}}^{4i*}Z_{\tilde{E}}^{4j}) \nonumber \\&-\upsilon \cos \beta |Y_{e_4}|^2\delta _{ij}-\frac{A_{E_4}}{\sqrt{2}}Z_{\tilde{E}}^{2i*}Z_{\tilde{E}}^{1j}\nonumber \\&-\frac{\mu ^*}{\sqrt{2}}Y_{e_5}Z_{\tilde{E}}^{4i*}Z_{\tilde{E}}^{3j}-\frac{1}{2}Y_{e_4}^*Z_{\tilde{E}}^{2j}\lambda _LZ_{\tilde{E}}^{4i*}\bar{\upsilon }_{NL} \nonumber \\&+\frac{1}{2}Z_{\tilde{E}}^{1i*}Y_{e_4}^*\lambda _EZ_{\tilde{E}}^{3j}\upsilon _{NL} \Big )\sin \alpha \Big ]. \end{aligned}$$
(16)

We also deduce the Z-exotic slepton–exotic slepton (\(Z-\tilde{E}-\tilde{E}\)) coupling in the EBLMSSM, which is given out as

$$\begin{aligned} \mathcal {L}_{Z\tilde{E}\tilde{E}}= & {} \sum _{i,j=1}^4\tilde{E}^{i*}\tilde{E}^{j}Z\frac{e}{2s_Wc_W} \nonumber \\&\quad \times \Big [\Big (Z_{\tilde{E}}^{1i*}Z_{\tilde{E}}^{1j}+Z_{\tilde{E}}^{4i*}Z_{\tilde{E}}^{4j}\Big ) -2s_W^2\delta _{ij}\Big ]. \end{aligned}$$
(17)

As the new introduced superfield in the EBLMSSM, Y leads to new couplings. The lepton–exotic lepton-Y coupling used in this work is shown here

$$\begin{aligned} \mathcal {L}_{lL'Y}= & {} \sum _{i,j=1}^2\bar{l}^I\Big (\lambda _4W_L^{1i}Z_Y^{1j*}P_R -\lambda _6U_L^{2i}Z_Y^{2j*}P_L\Big )\nonumber \\&\times L'_{i+3}Y_j^*+h.c. \end{aligned}$$
(18)

Superfield \(\tilde{Y}\) is also a new term beyond the BLMSSM. We deduce the lepton–exotic slepton-\(\tilde{Y}\) coupling as

$$\begin{aligned} \mathcal {L}_{l\tilde{E}\tilde{Y}}=\bar{\tilde{Y}}\Big (\lambda _4Z_{\tilde{E}}^{4i*}P_L -\lambda _6Z_{\tilde{E}}^{3i*}P_R\Big )l^I\tilde{E}_i^*+h.c. \end{aligned}$$
(19)

In the EBLMSSM, the new effects are added from the couplings of lepton–slepton–lepton neutralino, \(h^0(Z)\)-slepton–slepton, \(h^0(Z)\)-sneutrino–sneutrino, \(h^0(Z)\)-exotic lepton–exotic lepton and \(h^0(Z)\)-exotic neutrino–exotic neutrino, which are different from those in the BLMSSM. However, these couplings possess the same writing forms as those in the BLMSSM.

3 The processes \(l_j\rightarrow l_i \gamma \), muon conversion to electron in nuclei, the \(\tau \) decays and \(h^0\rightarrow l_i l_j\) in the EBLMSSM

In this section, we analyze the branching ratios of CLFV processes \(l_j\rightarrow l_i \gamma \), muon conversion rates to electron in Au nuclei, the branching ratios of rare \(\tau \) decays and \(h^0\rightarrow l_i l_j\) in the EBLMSSM.

3.1 Rare decays \(l_j\rightarrow l_i \gamma \)

Generally, the corresponding effective amplitude for processes \(l_j\rightarrow l_i\gamma \) can be written as [46]

$$\begin{aligned}&\mathcal{M}=e\epsilon ^{\mu }{\bar{u}}_i(p+q)[q^2\gamma _{\mu }(C_1^LP_L+C_1^RP_R) \nonumber \\&\quad \quad \quad +m_{l_j}i\sigma _{\mu \nu }q^{\nu }(C_2^LP_L+C_2^RP_R)]u_j(p), \nonumber \\&C_{\alpha }^{L,R}=C_{\alpha }^{L,R}(n)+C_{\alpha }^{L,R}(c)+C_{\alpha }^{L,R}(W),\alpha =1,2, \end{aligned}$$
(20)

where p (q) represents the injecting lepton (photon) momentum. \(m_{l_j}\) is the jth generation lepton mass. \(\epsilon \) is the photon polarization vector and \(u_i(p)\) (\(v_i(p)\)) is the lepton (antilepton) wave function. In Fig. 1, we show the relevant Feynman diagrams corresponding to above amplitude. The Wilson coefficients \(C_{\alpha }^{L,R}(\alpha =1,2)\) are discussed as follows.

\(C_{\alpha }^{L,R}(n)(\alpha =1,2)\), the virtual neutral fermion contributions corresponding to Fig. 1a, are deduced in the following form,

$$\begin{aligned} C_1^L(n)= & {} \sum _{F=\chi ^0/\chi _L^0,\nu ,\tilde{Y}}\sum _{S=\tilde{L},H^{\pm },\tilde{E}} \frac{1}{6m_{\Lambda }^2} H_R^{SF\bar{l}_i}H_L^{S^*l_j\bar{F}}I_4(x_F,x_S), \nonumber \\ C_2^L(n)= & {} \sum _{F=\chi ^0/\chi _L^0,\nu ,\tilde{Y}}\sum _{S=\tilde{L}, H^{\pm },\tilde{E}} \nonumber \\&\times \frac{m_F}{m_{l_j}m_{\Lambda }^2} H_L^{SF\bar{l}_i}H_L^{S^*l_j\bar{F}}[I_2(x_F,x_S)-I_3(x_F,x_S)], \nonumber \\ C_{\alpha }^R(n)= & {} C_{\alpha }^L(n)|_{L{\leftrightarrow }R},\alpha =1,2, \end{aligned}$$
(21)

where \(x_i=m_i^2/m_{\Lambda }^2\), \(m_i\) is the corresponding particle mass and \(m_{\Lambda }\) is the new physics energy scale. \(H_{L,R}^{SF\bar{l}_i}\) represent the left (right)-hand part of the coupling vertex. The concrete expressions for one-loop functions \(I_i(x_1,x_2)(i=1,2,3,4,5)\) are collected in “Appendix”.

Then, we discuss the virtual charged fermion contributions \(C_{\alpha }^{L,R}(c)(\alpha =1,2)\) corresponding to Fig. 1b

$$\begin{aligned} C_1^L(c)= & {} \sum _{F=\chi ^{\pm },L'}\sum _{S=\tilde{\nu },Y}\frac{1}{6m_{\Lambda }^2} H_R^{SF\bar{l}_i}H_L^{S^*l_j\bar{F}} \nonumber \\&\times [3I_2(x_S,x_F)-I_4(x_S,x_F)], \nonumber \\ C_2^L(c)= & {} \sum _{F=\chi ^{\pm },L'}\sum _{S=\tilde{\nu },Y}\frac{m_F}{m_{l_j}m_{\Lambda }^2} H_L^{SF\bar{l}_i}H_L^{S^*l_j\bar{F}}I_2(x_S,x_F), \nonumber \\ C_{\alpha }^R(c)= & {} C_{\alpha }^L(c)|_{L{\leftrightarrow }R},\alpha =1,2 \end{aligned}$$
(22)

Furthermore, the corrections from Fig. 1c are denoted by \(C_{\alpha }^{L,R}(W)(\alpha =1,2)\)

$$\begin{aligned} C_1^L(W)= & {} \sum _{F=\nu }\frac{1}{m_{\Lambda }^2} H_L^{WF\bar{l}_i}H_L^{W^*l_j\bar{F}} \nonumber \\&\times \left[ -2I_2(x_F,x_W)+\frac{1}{3}I_4(x_F,x_W)\right] , \nonumber \\ C_2^L(W)= & {} \sum _{F=\nu }\frac{1}{m_{\Lambda }^2} H_L^{WF\bar{l}_i}H_L^{W^*l_j\bar{F}}\frac{m_{l_i}}{m_{l_j}} \nonumber \\&\times [I_2(x_F,x_W)+I_4(x_F,x_W)], \nonumber \\ C_1^R(W)= & {} 0, \nonumber \\ C_2^R(W)= & {} \sum _{F=\nu }\frac{1}{m_{\Lambda }^2} H_L^{WF\bar{l}_i}H_L^{W^*l_j\bar{F}} \nonumber \\&\times [I_2(x_F,x_W)+I_4(x_F,x_W)]. \end{aligned}$$
(23)

However, the contributions from WW-neutrino diagram can be ignored due to the tiny neutrino mass.

We deduce the decay widths for processes \(l_j\rightarrow l_i \gamma \)

$$\begin{aligned} \Gamma \left( l_j\rightarrow l_i\gamma \right) =\frac{e^2}{16\pi }m_{l_j}^{5}\left( |C_2^L|^2+|C_2^R|^2\right) . \end{aligned}$$
(24)

Then, the concrete branching ratios of \(l_j\rightarrow l_i \gamma \) can be expressed as [46]

$$\begin{aligned} Br\left( l_j\rightarrow l_i\gamma \right) =\Gamma \left( l_j\rightarrow l_i\gamma \right) /\Gamma _{l_j}. \end{aligned}$$
(25)

Here, \(\Gamma _{l_j}\) represent the total decay widths of the charged leptons \(l_j\). We take \(\Gamma _{\mu }\simeq 2.996\times 10^{-19}\) GeV and \(\Gamma _{\tau }\simeq 2.265\times 10^{-12}\) GeV [18] in our latter numerical calculations.

Fig. 2
figure 2

The penguin type diagrams for the \(\mu -e\) conversion processes at the quark level

Fig. 3
figure 3

The box type diagrams for the \(\mu -e\) conversion processes at the quark level

3.2 \(\mu -e\) conversion in Au nuclei within the EBLMSSM

In this section, we just give out the figures for \(\mu -e\) conversion in nuclei at the quark level within the EBLMSSM, which are shown in Figs. 2 and 3. In the BLMSSM, the theoretical results for muon conversion to electron rates in nuclei are discussed specifically in our previous work [30, 33]. We find that Au nuclei currently give the most stringent bound on conversion rates, so we only study the \(\mu -e\) conversion rates in Au nuclei in this work. The new corrected particles in the EBLMSSM play important roles to this \(\mu -e\) conversion processes. Considering the constraints from \(\mu \rightarrow e\gamma \) within EBLMSSM, we study \(\mu -e\) conversion in Au nuclei, and the corresponding numerical results will be discussed in B of Sect. 4.

3.3 Rare \(\tau \) decays within the EBLMSSM

In this section, we discuss the rare \(\tau \) decays, which are \(\tau \rightarrow 3l_i\) and \(l_i\) represents particle e or \(\mu \). We give out both the penguin type diagrams and box type diagrams in Figs. 4 and 5. The theoretical results for the \(\tau \) decays are discussed specifically in our previous work [32]. In the EBLMSSM, the numerical results of \(\tau \) decays can be influenced by the new corrected particles, such as exotic lepton (slepton), slepton (sneutrino), lepton neutralino, Y and \(\tilde{Y}\). In our latter work, we will analyze this \(\tau \) decays in detail.

Fig. 4
figure 4

The penguin type diagrams for the \(\tau \) decays

Fig. 5
figure 5

The box type diagrams for the \(\tau \) decays

3.4 Rare decay \(h^0\rightarrow l_i l_j\)

The corresponding effective amplitude for \(h^0\rightarrow \bar{l}_i l_j\) can be summarized as

$$\begin{aligned}&\mathcal{A}={\bar{u}}_i(q)(N_LP_L+N_RP_R)v_j(p), \nonumber \\&N_{L,R}=N_{L,R}(S_1)+N_{L,R}(S_2)+N_{L,R}(W) \nonumber \\&\quad \quad \quad +A_{L,R}(S_1)+A_{L,R}(S_2)+A_{L,R}(W_1)+A_{L,R}(W_2).\nonumber \\ \end{aligned}$$
(26)

Here \(N_{L,R}(S_1)\) are the coupling coefficients corresponding to triangle diagrams in Fig. 6a, \(N_{L,R}(S_2)\) denote the contributions from Fig. 6b. The effects from Fig. 6c, d can be shown by \(N_{L,R}(W)\). \(A_{L,R}(S_1)\) and \(A_{L,R}(S_2)\) represent the contributions from self-energy diagrams Figs.7a, and 6b respectively. The effects from Fig. 7c, d can be summarized by \(A_{L,R}(W_1)\) and \(A_{L,R}(W_2)\) respectively. We give out the concrete expressions for these contributions as follows.

Fig. 6
figure 6

The triangle type diagrams for decays \(h^0\rightarrow \bar{l}_i l_j\)

Fig. 7
figure 7

The self-energy type diagrams for decays \(h^0\rightarrow \bar{l}_i l_j\)

The contributions from triangle diagrams in Fig. 6:

$$\begin{aligned} N_L(S_1)= & {} \sum _{F=\chi ^{\pm }\chi ^0/\chi _L^0,\nu ,\tilde{Y}}\sum _{S=\tilde{\nu },\tilde{L},H^{\pm },\tilde{E}} \nonumber \\&\times \frac{m_F}{m_{\Lambda }^2}H_L^{S_2F\bar{l}_i}H^{h^0S_1S_2^*}H_L^{S_1^*l_j\bar{F}}G_1(x_F,x_{S_1},x_{S_2}), \nonumber \\ N_R(S_1)= & {} N_L(S_1)|_{L{\leftrightarrow }R}, \end{aligned}$$
(27)
$$\begin{aligned} N_L(S_2)= & {} \sum _{F=\chi ^{\pm },\chi ^0,\nu ,L'}\sum _{S=\tilde{\nu },\tilde{L},H^{\pm },Y} \nonumber \\&\times \Bigg [H_L^{SF_2\bar{l}_i}H_R^{h^0F_1\bar{F}_2}H_L^{S^*l_j\bar{F}_1}G_2(x_S,x_{F_1},x_{F_2}) \nonumber \\&+\,\frac{m_{F_1}m_{F_2}}{m_{\Lambda }^2} H_L^{SF_2\bar{l}_i}H_L^{h^0F_1\bar{F}_2}H_L^{S^*l_j\bar{F}_1}G_1(x_S,x_{F_1},x_{F_2})\Bigg ], \nonumber \\ N_R(S_2)= & {} N_L(S_2)|_{L{\leftrightarrow }R}, \end{aligned}$$
(28)
$$\begin{aligned} N_L(W)= & {} -\sum _{F=\nu }\frac{m_{l_i}}{m_{\Lambda }^2}H_L^{WF\bar{l}_i}H_L^{\bar{F}l_jW^*}H^{h^0WW^*}I_2(x_F,x_W) \nonumber \\&+\,\sum _{F_1,F_2=\nu }\sqrt{x_{l_i}x_{F_2}}H_L^{WF_2\bar{l}_i}H_L^{h^0F_1\bar{F}_2}H_L^{\bar{F}_1l_jW^*} \nonumber \\&\times [G_1(x_W,x_{F_1},x_{F_2})+x_{F_2}\frac{d}{dx_{F_2}}G_1(x_W,x_{F_1},x_{F_2})], \nonumber \\ N_R(W)= & {} -\sum _{F=\nu }\frac{m_{l_j}}{m_{\Lambda }^2}H_L^{WF\bar{l}_i}H_L^{\bar{F}l_jW^*}H^{h^0WW^*}I_2(x_F,x_W) \nonumber \\&+\,\sum _{F_1,F_2=\nu }\sqrt{x_{l_j}x_{F_2}}H_L^{WF_2\bar{l}_i}H_L^{h^0F_1\bar{F}_2}H_L^{\bar{F}_1l_jW^*} \nonumber \\&\times \left[ G_1(x_W,x_{F_1},x_{F_2})+x_{F_1}\frac{d}{dx_{F_1}}G_1(x_W,x_{F_1},x_{F_2})\right] .\nonumber \\ \end{aligned}$$
(29)

The contributions from self-energy type diagrams correspond to Fig. 7:

$$\begin{aligned}&A_L(S_1)=\sum _{F=\chi ^0/\chi _L^0,\chi ^{\pm },L',\tilde{Y}}\sum _{S=\tilde{L},\tilde{\nu },Y,\tilde{E}} \nonumber \\&\quad \times \frac{1}{m_{l_j}^2-m_{l_i}^2} H^{h^0l_i\bar{l}_i}\left\{ m_F\left( m_{l_j}H_R^{\bar{l}_iFS}H_R^{l_j\bar{F}S^*} \right. \right. \nonumber \\&\quad \left. +\,m_{l_i}H_L^{\bar{l}_iFS}H_L^{l_j\bar{F}S^*}\right) \nonumber \\&\quad \times \left[ I_1(x_F,x_S)+\frac{m_{l_j}^2}{m_{\Lambda }^2}(I_2(x_F,x_S)-I_3(x_F,x_S))\right] \nonumber \\&\quad -\frac{1}{2}\left( m_{l_j}^2H_R^{\bar{l}_iFS}H_L^{l_j\bar{F}S^*} +m_{l_i}m_{l_j}H_L^{\bar{l}_iFS}H_R^{l_j\bar{F}S^*}\right) \nonumber \\&\quad \left. I_5(x_F,x_S)\phantom {\int }\right\} , \nonumber \\&A_R(S_1)=A_L(S_1)|_{L{\leftrightarrow }R}, \end{aligned}$$
(30)
$$\begin{aligned}&A_L(S_2)=\sum _{F=\chi ^0/\chi _L^0,\chi ^{\pm },L',\tilde{Y}}\sum _{S=\tilde{L},\tilde{\nu },Y,\tilde{E}} \nonumber \\&\quad \frac{1}{m_{l_i}^2-m_{l_j}^2} H^{h^0l_j\bar{l}_j}\left\{ \phantom {\left. -\frac{1}{2}(m_{l_i}^2H_L^{\bar{l}_iFS}H_R^{l_j\bar{F}S^*}+m_{l_i}m_{l_j}H_R^{\bar{l}_iFS}H_L^{l_j\bar{F}S^*}) I_5(x_F,x_S)\right\} }m_F\left( m_{l_i}H_R^{\bar{l}_iFS}H_R^{l_j\bar{F}S^*} +m_{l_j}H_L^{\bar{l}_iFS}H_L^{l_j\bar{F}S^*}\right) \right. \nonumber \\&\quad \times \left[ I_1(x_F,x_S)+\frac{m_{l_i}^2}{m_{\Lambda }^2}(I_2(x_F,x_S)-I_3(x_F,x_S))\right] \nonumber \\&\quad \left. -\frac{1}{2}(m_{l_i}^2H_L^{\bar{l}_iFS}H_R^{l_j\bar{F}S^*}+m_{l_i}m_{l_j}H_R^{\bar{l}_iFS}H_L^{l_j\bar{F}S^*}) I_5(x_F,x_S)\right\} , \nonumber \\&A_R(S_2)=A_L(S_2)|_{L{\leftrightarrow }R}, \end{aligned}$$
(31)
$$\begin{aligned}&A_L(W_1)=-\sum _{F=\nu }\frac{m_{l_j}^2}{m_{l_j}^2-m_{l_i}^2} H^{h^0l_i\bar{l}_i}H_L^{\bar{l}_iFW}H_L^{l_j\bar{F}W^*}I_5(x_F,x_W), \nonumber \\&A_R(W_1)=-\sum _{F=\nu }\frac{m_{l_i}m_{l_i}}{m_{l_j}^2-m_{l_i}^2} H^{h^0l_j\bar{l}_j}H_L^{\bar{l}_iFW}H_L^{l_j\bar{F}W^*}I_5(x_F,x_W). \end{aligned}$$
(32)
$$\begin{aligned}&A_L(W_2)=-\sum _{F=\nu }\frac{m_{l_i}m_{l_j}}{m_{l_i}^2-m_{l_j}^2} H^{h^0l_j\bar{l}_j}H_L^{\bar{l}_iFW}H_L^{l_j\bar{F}W^*}I_5(x_F,x_W), \nonumber \\&A_R(W_2)=-\sum _{F=\nu }\frac{m_{l_i}^2}{m_{l_i}^2-m_{l_j}^2} H^{h^0l_j \bar{l}_j}H_L^{\bar{l}_iFW}H_L^{l_j\bar{F}W^*}I_5(x_F,x_W), \end{aligned}$$
(33)

where, the one-loop functions \(G_i(x_1,x_2,x_3)(i=1,2)\) are collected in “Appendix”.

The decay widths for processes \(h^0\rightarrow l_i l_j\) are deduced here

$$\begin{aligned} \Gamma (h^0\rightarrow l_il_j)=(h^0\rightarrow \bar{l}_il_j)+(h^0\rightarrow l_i\bar{l}_j), \end{aligned}$$
(34)

where \(\Gamma \left( h^0\rightarrow \bar{l}_il_j\right) =\frac{1}{16\pi }m_{h^0}\left( |N_L|^2+|N_R|^2\right) \)[47, 48]. Correspondingly, the calculations for \(h^0\rightarrow l_i\bar{l}_j\) are same as those for \(h^0\rightarrow \bar{l}_il_j\).

Above all, the branching ratios of \(h^0\rightarrow l_il_j\) can be summarized as

$$\begin{aligned} Br(h^0\rightarrow l_il_j)=\Gamma (h^0\rightarrow l_il_j)/\Gamma _{h^0}. \end{aligned}$$
(35)

Here, the total decay width of the 125.1 GeV Higgs boson is \(\Gamma _{h^0}\simeq 4.1\times 10^{-3}\) GeV [18].

\(B^0\) meson is made up of d \(\bar{b}\) and \(B_s^0\) meson is constituted of s \(\bar{b}\). The present experiment upper bounds for \(B^0\) and \( B_s^0\) meson decays are respectively \(Br(B^0\rightarrow e^+\mu ^-)<2.8\times 10^{-9}\) and \(Br(B_s^0\rightarrow e^+\mu ^-)<1.1\times 10^{-8}\) [18]. New contributions to rare \(B^0\) and \( B_s^0\) meson decays emerge at one-loop level with the box diagrams. In the EBLMSSM, the redefined particles sleptons and sneutrinos lead to new effects to these rare \(B^0\) and \( B_s^0\) meson decays. So parameters \(\tan \beta _{NL}\) and \(v_{Nlt}\) may play the dominated roles to the \(B^0\) and \(B_s^0\) meson decays. \(\pi ^+,K^+\) mesons are respectively comprised of u \(\bar{d}\) and u \(\bar{s}\). Particle Date Group gives us the present experiment upper bounds for \((\pi ^+/K^+)\rightarrow l_i^+\nu _j\), which are \(Br(\pi ^+\rightarrow \mu ^+\nu _e)<8.0\times 10^{-3}\) and \(Br(K^+\rightarrow \mu ^+\nu _e)<4.0\times 10^{-3}\) [18]. In the EBLMSSM, the penguin type diagrams, self-energy type diagrams and box type diagrams all affect the processes \((\pi ^+/K^+)\rightarrow l_i^+\nu _j, i\ne j\). CLFV contributions arise from loop corrections with the \(W^{\pm }\) and heavy charged Higgs propagator. Furthermore, the loop contributions are also related with the exotic slepton (sneutrino), exotic lepton (neutrino), lepton neutralino and slepton (sneutrino) particles. Therefore, processes \((\pi ^+/K^+)\rightarrow l_i^+\nu _j\) will be strongly affected by parameters presented in the EBLMSSM. We hope a detailed analysis is going to be discussed in our next work.

4 Numerical results

In this section, we discuss the numerical results. In our previous work [34], we research the processes \(h^0\rightarrow \gamma \gamma \), \(h^0\rightarrow VV, V=(Z,W)\) in the EBLMSSM, and the corresponding numerical results are discussed in Sect. 5.1 of work [34]. The CP-even Higgs masses \(m_{h^0}, m_{H^0}\) and CP-odd Higgs mass \(m_A^0\) are also analyzed. In the reasonable parameter space, the values of branching ratios for \(h^0\rightarrow \gamma \gamma \)(\(R_{\gamma \gamma }\)) and \(h^0\rightarrow VV\)(\(R_{VV}\)) both meet the experiment limits. Therefore, the Higgs decays in the EBLMSSM play important roles to promote physicists to explore new physics. And the corresponding constraints are also considered in our work. The CP-even Higgs mass is considered as an input parameter, which is \(m_{h^0}=125.1\) GeV in our latter numerical discussions.

In the EBLMSSM, to obtain a more transparent numerical results, we adopt the following assumptions on parameter space:

$$\begin{aligned}&Y_{u_4} = 1.2Y_t,~Y_{u_5} = 0.6Y_t,~ Y_{d_4}=Y_{d_5} = 2Y_b, \nonumber \\&\quad Y_{\nu _4} = Y_{\nu _5} = 0.8,~\mu _B = \mu _L = 0.5 \mathrm{TeV},\nonumber \\&m_{\tilde{Q}_4} = m_{\tilde{Q}_5} =m_{\tilde{U}_4} = m_{\tilde{U}_5} = m_{\tilde{D}_4} \nonumber \\&\quad \quad \quad = m_{\tilde{D}_5} = m_{\tilde{\nu }_4} = m_{\tilde{\nu }_5} =1\mathrm{TeV},~B_4 = L_4 = 1.5,\nonumber \\&A_{u_4} =A_{u_5} = A_{d_4} = A_{d_5} =A_{\nu _4} = A_{\nu _5} =1\mathrm{TeV}, \nonumber \\&\quad \lambda _Q = 0.4,\lambda _u=\lambda _d = 0.5,(\lambda _{Nc})_{ii} = 1,\nonumber \\&A_{BQ}=A_{BU} = A_{BD} = 1\mathrm{TeV},~g_B = 1/3,~g_L = 1/6, \nonumber \\&\quad \tan \beta _B =1.5,~\tan \beta _L = 2,\nonumber \\&m_{\tilde{Q}_3}=m_{\tilde{U}_3}= 1.2\mathrm{TeV}, \nonumber \\&\quad m_{\tilde{D}_3}= 1.5\mathrm{TeV},~A_t = 1.7\mathrm{TeV}, \nonumber \\&\quad A_b = 3\mathrm{TeV},~M_L = 1\mathrm{TeV}, \nonumber \\&(m_{\tilde{\nu }})_{ii}= 1\mathrm{TeV}, \nonumber \\&\quad (A_{N})_{ii} =(A_{Nc})_{ii} = 0.5\mathrm{TeV},~~m_{\Lambda }=1\mathrm{TeV} \end{aligned}$$
(36)

where \(i=1,2,3\), \(Y_t\) (\(Y_b\)) corresponds to the Yukawa coupling constant of top (bottom) quark, whose concrete form can be written as \(Y_t = \sqrt{2} m_t/(\upsilon \sin \beta )\) (\(Y_b = \sqrt{2} m_b/(\upsilon \cos \beta )\)).

In order to simplify the numerical analysis, we use the following assumptions:

$$\begin{aligned}&m_{\tilde{L}_4} = m_{\tilde{L}_5} =m_{\tilde{e}_4} = m_{\tilde{e}_5} = M_{\tilde{E}},\nonumber \\&A_{e_4} =A_{e_5} = A_{\tilde{E}},\lambda _L =\lambda _E=\lambda _{N_L} = L_l,\nonumber \\&A_{LL}=A_{LE} = A_{LN} = A_E,(\lambda _4^2)^{IJ} \nonumber \\&\quad =(\lambda _6^2)^{IJ} = (Lm^2)^{IJ},I,J=1,2,3,v_{Nlt}=v_N,\nonumber \\&(m_{\tilde{L}}^2)_{ii}=(M_{Ls})_{ii}^2,~(m_{\tilde{L}}^2)_{ij} =M_{Lf}^2,~(A_l)_{ii}=Al, \nonumber \\&\quad (A'_l)_{ii}=A'l,i,j=1.2.3,i\ne j. \end{aligned}$$
(37)

We take \(\sqrt{(Lm^2)^{12}}=L_F\) and \(\sqrt{(Lm^2)^{13}}=\sqrt{(Lm^2)^{23}}=L_f\).

4.1 \(l_j\rightarrow l_i \gamma \)

4.1.1 \(\mu \rightarrow e \gamma \)

CLFV process \(\mu \rightarrow e \gamma \) contributes to explore the new physics, whose experiment upper bound of the branching ratio is around \(5.7 \times 10^{-13}\) at \(90\%\) confidence level. In this part, we discuss the effects on process \(\mu \rightarrow e \gamma \) from some new introduced parameters in the EBLMSSM.

Parameter \(A_{\tilde{E}}\) is present in the non-diagonal parts of the exotic slepton mass matrix. Parameter \(Y_{e5}\) is not only related to the non-diagonal parts of the exotic lepton and exotic slepton mass matrices, but also connected with the diagonal exotic slepton elements. In the EBLMSSM, exotic lepton and exotic slepton are both different from those in the BLMSSM. We assume that \(M_{\tilde{E}}=\mu _Y=1.5\) TeV, \(A_E=\mu _{NL}=1\) TeV, \(L_l=1\), \(L_F=10^{-3}\), \(\mu =0.7\) TeV, \(B_Y=0.94\) TeV, \(Y_{e4}=0.5\), \((M_{Ls})_{11}^2=6 \mathrm{TeV}^2\), \((M_{Ls})_{22}^2=4 \mathrm{TeV}^2\), \((M_{Ls})_{33}^2=1 \mathrm{TeV}^2\), \(M_{Lf}^2=10^{-3} \mathrm{TeV}^2\), \(Al=2\) TeV, \(A'l=0.3\) TeV, \(m_1=m_2=1.5\mathrm{TeV}\), \(\tan \beta =6\) and \(\tan \beta _{NL}=2\) . With \(Y_{e5}=1.0(1.5,2.0)\), the branching ratios of \(\mu \rightarrow e \gamma \) versus parameter \(A_{\tilde{E}}\) are studied, which are shown in Fig. 8. When \(A_{\tilde{E}}\) is in the region \(0.1 \sim 2.5\) TeV, the numerical results change from \(5\times 10^{-15}\) to \(5\times 10^{-13}\). These three lines all increase quickly and approach to the experiment upper bound. Therefore, \(A_{\tilde{E}}\) affects the numerical results strongly. Furthermore, corresponding to same \(A_{\tilde{E}}\), the solid line results are about two times as the dashed line results, and the dashed line results are about two times as the dotted line results. Larger \(Y_{e5}\) can lead to larger numerical results.

Fig. 8
figure 8

With \(Y_{e5}=1.0(1.5,2.0)\), the branching ratios of \(\mu \rightarrow e \gamma \) versus parameter \(A_{\tilde{E}}\) are plotted by the dotted, dashed and solid lines respectively

Fig. 9
figure 9

With \(B_Y=0.4(0.8,1.2)\) TeV, the branching ratios of \(\mu \rightarrow e \gamma \) versus parameter \(\mu _Y\) are plotted by the dotted, dashed and solid lines respectively

With the introduced superfields Y and \(Y'\) in the EBLMSSM, we deduce the Y and \(\tilde{Y}\) mass matrices. Parameters \(\mu _Y\) and \(B_Y\) are respectively present in the diagonal and non-diagonal terms of the Y mass matrix. And the mass of \(\tilde{Y}\) possesses the same value as \(\mu _Y\). So these two parameters affect the Y-lepton-exotic lepton and \(\tilde{Y}\)-lepton–exotic slepton couplings. Furthermore, these new couplings make contributions to the numerical results. Using \(Y_{e4}=0.8\), \(Y_{e5}=1.5\) and \(A_{\tilde{E}}=1\) TeV, we plot the branching ratios changing with \(\mu _Y\) in Fig. 9. The dotted (dashed, solid) line represents \(B_Y=0.4(0.8,1.2)\) TeV. We find that the branching ratios decrease quickly with the increasing \(\mu _Y\), which indicates that the large \(\mu _Y\) can restrain the numerical results evidently. Furthermore, the numerical results of these three lines are almost same with the unchanging \(\mu _Y\), so the contributions from parameter \(B_Y\) is small.

Fig. 10
figure 10

With \(\mu _{NL}=0.7(1.0,1.3)\) TeV, the branching ratios of \(\mu \rightarrow e \gamma \) versus parameter \(A_E\) are plotted by the dotted, dashed and solid lines respectively

Then, we study effects from the parameters \(A_E\) and \(\mu _{NL}\) on our numerical results. In EBLMSSM, parameters \(A_E\) and \(\mu _{NL}\) are both the non-diagonal elements in the exotic slepton and exotic sneutrino mass matrices. \(\mu _{NL}\) is also the non-diagonal element of lepton neutralino mass matrix. In Fig. 10, we present the branching ratios of \(\mu \rightarrow e \gamma \) versus \(A_E\) with \(\mu _{NL}=0.7(1.0,1.3)\)TeV, and the concrete results are plotted by dotted (dashed, solid) line. These three lines all increase quickly when \(A_E\) ranges from 0.1 to 1.8 TeV. Therefore, as the sensitive parameters in the EBLMSSM, the large \(A_E\) produces the large contributions on the results. However, the numerical results slightly decrease with the enlarging \(\mu _{NL}\), and the effects from \(\mu _{NL}\) are not so obvious as that \(A_E\).

4.1.2 \(\tau \rightarrow \mu \gamma \) (\(\tau \rightarrow e\gamma \))

In a similar way, the CLFV processes \(\tau \rightarrow e \gamma \) and \(\tau \rightarrow \mu \gamma \) are studied. The corresponding experimental upper bounds of the branching ratios are \(Br(\tau \rightarrow e \gamma )<3.3\times 10^{-8}\) and \(Br(\tau \rightarrow \mu \gamma )<4.4\times 10^{-8}\).

As a new introduced parameter in the EBLMSSM, parameter \(v_{N}\) is present in the mass matrices of slepton, sneutrino, exotic lepton (neutrino), exotic slepton (sneutrino) and lepton neutralino. In this part, we research the branching ratios of \(\tau \rightarrow e \gamma \) and \(\tau \rightarrow \mu \gamma \) changing with \(v_{N}\). Supposing \(M_{\tilde{E}}=\mu _Y=1.5\) TeV, \(A_E=\mu _{NL}=1\) TeV, \(A_{\tilde{E}}=1\) TeV, \(L_l=1\), \(L_f=0.08\), \(\mu =0.7\) TeV, \(B_Y=0.94\) TeV, \(Y_{e4}=Y_{e5}=0.8\), \((M_{Ls})_{ii}^2=S_m^2=1 \mathrm{TeV}^2\), i \(=\) 1,2,3, \(M_{Lf}^2=10^{-2} \mathrm{TeV}^2\), \(Al=2\) TeV, \(A'l=0.3\) TeV, \(m_1=m_2=1.5\mathrm{TeV}\), \(\tan \beta =6\) and \(\tan \beta _{NL}=2\), we plot the numerical results with \(v_{N}\) in Fig. 11 by dashed line and solid line respectively. Obviously, when the values of \(v_{N}\) change from 1.5 to 3.5 TeV, the results of \(Br(\tau \rightarrow e \gamma )\) and \(Br(\tau \rightarrow \mu \gamma )\) both shrink quickly. This implies that \(v_{N}\) is a sensitive parameter. Though the figure of process \(\tau \rightarrow e \gamma \) is under that of \(\tau \rightarrow \mu \gamma \), the both lines possess almost the same results when \(v_{N}\) takes same value. So we only study the branching ratios of process \(\tau \rightarrow \mu \gamma \) in following discussion.

Fig. 11
figure 11

The branching ratios of \(\tau \rightarrow e \gamma \) and \(\tau \rightarrow \mu \gamma \) versus parameter \(v_{N}\) are plotted by the dashed line and solid line respectively

Fig. 12
figure 12

With \(L_l=0.7(1.0,1.3)\), the branching ratios of \(\tau \rightarrow \mu \gamma \) versus parameter \(M_{\tilde{E}}^2\) are plotted by the dotted, dashed and solid lines respectively

Appearing in the diagonal terms of the exotic slepton and exotic sneutrino mass squared matrices, \(M_{\tilde{E}}\) affects the \(\tilde{Y}\)-lepton–exotic slepton coupling. Parameter \(L_l\), not only in the exotic lepton (neutrino) but also in exotic slepton (sneutrino) mass matrices, produces contributions to the numerical results through Y-lepton–exotic lepton and \(\tilde{Y}\)-lepton-exotic slepton couplings. With \((M_{Ls})_{11}^2=6 \mathrm{TeV}^2\), \((M_{Ls})_{22}^2=4 \mathrm{TeV}^2\), \((M_{Ls})_{33}^2=1 \mathrm{TeV}^2\) and \(v_{N}=3\) TeV, the numerical results versus \(M_{\tilde{E}}^2\) are plotted in Fig. 12. The dotted (dashed, solid) line corresponds to \(L_l=0.7(1.0,1.3)\). The figure shows that these three lines all decrease quickly when \(M_{\tilde{E}}^2\) varies from \(1\times 10^{6}\) to \(9\times 10^{6}\mathrm{GeV}^2\). With the same \(M_{\tilde{E}}^2\), the branching ratio decreases remarkably when \(L_l\) increases. Especially, the line is much steeper with \(L_l=0.7\) than that \(L_l=1.0(1.3)\). Obviously, both \(M_{\tilde{E}}\) and \(L_l\) are sensitive parameters to our numerical results.

Parameter \(L_f\) influences the numerical results through Y-lepton–exotic lepton and \(\tilde{Y}\)-slepton–exotic slepton couplings. And parameter \(Y_{e4}\) affects our numerical results through exotic lepton and exotic slepton. We discuss the numerical results with \(L_f\) varying from 0.01 to 0.3 in Fig. 13. The dotted (dashed, solid) line corresponds to \(Y_{e4}=0.5(1.0,1.5)\). The branching ratios possess slight changes when \(Y_{e4}\) takes different values for the unchanged \(L_f\), which indicates the effects from \(Y_{e4}\) can be ignored in our following discussion. It is easy to see that the numerical results increase sharply with the enlarging \(L_f\). So the non-diagonal elements of parameters \(\lambda _4^2\) and \(\lambda _6^2\) play important roles in our numerical studies.

Fig. 13
figure 13

With \(Y_{e4}=0.5(1.0,1.5)\), the branching ratios of \(\tau \rightarrow \mu \gamma \) versus parameter \(L_f\) are plotted by the dotted, dashed and solid lines respectively

4.2 \(\mu -e\) conversion rates in Au nuclei

The present sensitivity for the muon conversion rates to electron in Au nuclei is \(CR(\mu \rightarrow e:\;_{79}^{197}Au)<7\times 10^{-13}\). Considering the parameter constrains from \(\mu \rightarrow e \gamma \), we analyze the numerical results for this \(\mu -e\) conversion in Au nuclei.

Fig. 14
figure 14

With \(\tan \beta _{NL}=1.8,2.2,2.6\), the \(\mu -e\) conversion rates in Au nuclei versus parameter \(L_F\) are plotted by the dotted, dashed and solid lines respectively

Fig. 15
figure 15

With \(\mu =0.4(0.5,0.6)\) TeV, the \(\mu -e\) conversion rates in Au nuclei versus parameter \(M_{Lf}^2\) are plotted by the dotted, dashed and solid lines respectively

As the non-diagonal elements of matrix \((Lm^2)^{IJ}\) in the EBLMSSM, \(\sqrt{(Lm^2)^{12}}=L_F\) affects the numerical results through exotic lepton and exotic slepton. Choosing \(M_{\tilde{E}}=A_E=\mu _{NL}=\mu =1\) TeV, \(\mu _Y=2\) TeV, \(L_l=0.8\), \(B_Y=1.5\) TeV, \(Y_{e4}=0.8\), \(Y_{e5}=1.5\), \(S_m^2=6 \mathrm{TeV}^2\), \(M_{Lf}^2=500 \mathrm{GeV}^2\), \(m_1=m_2=3\mathrm{TeV}\), \(Al=2\) TeV, \(A'l=0.3\) TeV and \(\tan \beta =6\), we analyze the \(\mu -e\) conversion rates in Au nuclei with \(L_F\) in Fig. 14. \(\tan \beta _{NL}=1.8,2.2,2.6\) correspond to the dotted, dashed and solid lines respectively. When \(L_F\) changes from 0.001 to 0.006, These three lines all enlarge quickly and can easily reach the present sensitivity. So \(L_F\) greatly contributes to the numerical results. Furthermore, as \(L_F\) takes the same value, the larger \(\tan \beta _{NL}\), the smaller numerical result it is. The above analyses indicate \(L_F\) and \(\tan \beta _{NL}\) are both sensitive parameters.

As the non-diagonal elements of slepton and sneutrino mass matrices, \(M_{Lf}\) lead to strong mixing for slepton (sneutrino) with different generations. The parameter \(\mu \) presents in the mass matrices of slepton, sneutrino, exotic slepton and exotic sneutrino. So we study the \(\mu -e\) conversion rates in Au nuclei versus parameters \(M_{L_f}^2\). As \(M_{\tilde{E}}=2\) TeV, \(\mu _Y=2\) TeV, \(B_Y=0.9\) TeV, \(S_m^2=12 \mathrm{TeV}^2\), \(m_1=m_2=3\mathrm{TeV}\), \(Al=1.9\) TeV, \(\tan \beta _{NL}=2\) and \(L_F=0.001\), we show the numerical results changing with \(M_{Lf}^2\), which are given in Fig. 15. The dotted, dashed and solid lines respectively correspond to \(\mu =0.4,0.5,0.6\) TeV. With the enlarging \(M_{Lf}^2\), the numerical results increase quickly. As \(M_{Lf}^2>5000\,\mathrm{GeV}^2\) and taking the same values, these three lines almost possess similar results, which indicates that the effects from parameter \(\mu \) are small.

4.3 \(\tau \) decays

The experiment upper bounds for \(\tau \) decays are \(Br(\tau \rightarrow 3e)<2.7\times 10^{-8}\) and \(Br(\tau \rightarrow 3\mu )<2.1\times 10^{-8}\). Considering the constraints from \(\tau \rightarrow e \gamma \) and \(\tau \rightarrow \mu \gamma \), we discuss the numerical results for decays \(\tau \rightarrow 3e\) and \(\tau \rightarrow 3\mu \).

Using \(\mu _Y=3\) TeV, \(M_{\tilde{E}}=A_E=\mu _{NL}=1\) TeV, \(L_l=0.8\), \(\mu =0.7\) TeV, \(B_Y=1.5\) TeV, \(Y_{e4}=0.8\), \(Y_{e5}=1.5\), \((M_{Ls})_{ii}^2=6 \mathrm{TeV}^2\), \(M_{Lf}^2=500 \mathrm{GeV}^2\), \(Al=1\) TeV, \(A'l=0.3\) TeV, \(m_1=m_2=0.5\mathrm{TeV}\), \(\tan \beta =6\) and \(\tan \beta _{NL}=1.5\), we plot the numerical results of \(\tau \rightarrow 3e\) and \(\tau \rightarrow 3\mu \) in Fig. 16 by dotted line and solid line respectively. With parameter \(L_f\) changing from 0.01 to 0.3, the values of these two lines are almost the same and both increase quickly. So we just discuss the numerical results for \(\tau \rightarrow 3e\) decays as follows (Fig. 17).

Fig. 16
figure 16

The branching ratios of \(\tau \rightarrow 3\mu \) and \(\tau \rightarrow 3e\) versus parameter \(L_f\) are plotted by the solid line and dashed line respectively

Fig. 17
figure 17

With \(\tan \beta _{NL}=1.5(1.7,1.9)\), the branching ratios of \(\tau \rightarrow 3e\) versus parameter \(\mu _{NL}\) are plotted by the dotted, dashed and solid lines respectively

Choosing \(Y_{e4}=0.8\), \(Y_{e5}=1.5\), \(M_{\tilde{E}}=\mu _Y=1.5\) TeV and \(B_Y=0.9\) TeV, we study the branching ratios of \(\tau \rightarrow 3e\) changing with parameter \(\mu _{NL}\). The numerical results varying with \(\tan \beta _{NL}=1.5(1.7,1.9)\) are plotted by the dotted, dashed and solid lines respectively. As \(\mu _{NL}\) changes from 0.3 to 1 TeV, these three lines all have the obviously improvement. As \(\mu _{NL}>1\) TeV, the numerical results increase slowly. Besides, when \(\mu _{NL}\) dose not change, the branching ratios of \(\tau \rightarrow 3e\) enlarge with the increased \(\tan \beta _{NL}\), and the bigger \(\tan \beta _{NL}\), the bigger change it is in the graph. Therefore, both \(\mu _{NL}\) and \(\tan \beta _{NL}\) affect the numerical results in a certain degree.

4.4 \(h^0\rightarrow l_i l_j\)

In this part, we study the CLFV processes \(h^0\rightarrow l_i l_j\). The most strict constraint \(m_{h^0}=125.1\) GeV is considered as an input parameter. We also take into account the limits from processes \(l_j\rightarrow l_i \gamma \), the muon conversion to electron in Au nuclei and the \(\tau \) decays discussed above.

4.4.1 \(h^0\rightarrow \mu \tau (h^0\rightarrow e \tau )\)

At first, we picture the branching ratios of decays \(h^0\rightarrow \mu \tau \) and \(h^0\rightarrow e \tau \) versus \(A_{\tilde{E}}\) in Fig. 18. We choose the relevant parameters as \(\mu _Y=1.5\) TeV, \(M_{\tilde{E}}=1.4\) TeV, \(A_E=m_1=1\) TeV, \(\mu _{NL}=2\) TeV, \(L_l=1\), \((Lm^2)^{13}=(Lm^2)^{23}=L_f=0.3^2\), \(\mu =0.7\) TeV, \(B_Y=0.94\) TeV, \(Y_{e4}=0.8\), \(Y_{e5}=1.5\), \(S_m^2=1 \mathrm{TeV}^2\), \(M_{Lf}^2=12000 \mathrm{GeV}^2\), \(Al=1\) TeV, \(A'l=3\) TeV, \(m_2=0.5\) TeV, \(\tan \beta =6\) and \(\tan \beta _{NL}=2\). Although the line of \(h^0\rightarrow \mu \tau \) is under that of \(h^0\rightarrow e \tau \), these two processes almost have the same variation trend. With the enlarging \(A_{\tilde{E}}\), the numerical results increase quickly.

Fig. 18
figure 18

The branching ratios of \(h^0\rightarrow \mu \tau \) and \(h^0\rightarrow e \tau \) versus parameter \(A_{\tilde{E}}\) are plotted by the solid line and dashed line respectively

Fig. 19
figure 19

With \(S_m^2=1(2,3)\mathrm{TeV}^2\), the branching ratios of \(h^0\rightarrow \mu \tau \) versus parameter \(\tan {\beta }\) are plotted by the dotted, dashed and solid lines respectively

Then the effects from the parameters \(\tan {\beta }\) and \(S_m\) are studied. \(\tan {\beta }\) is related to \(v_u\) and \(v_d\), and appears in almost all mass matrices of CLFV processes. \(S_m\) are present in the diagonal elements of slepton and sneutrino mass matrices. With \(A_{\tilde{E}}=2\) TeV, \(L_f=0.25\), \(Y_{e4}=1.2\) and \(Y_{e5}=0.8\), Fig. 19 shows the branching fractions of \(h^0\rightarrow \mu \tau \) varying with the parameter \(\tan {\beta }\). \(S_m^2=1(2,3)\,\mathrm{TeV}^2\) corresponds to the dotted (dashed, solid) line. These three lines almost overlap, so the effects from \(S_m\) are small. As \(\tan {\beta }\) varies from 6 to 9, the numerical results decrease obviously. As \(\tan {\beta }>9\), the numerical results increase quickly. So \(\tan {\beta }\) plays very important roles to CLFV processes.

4.4.2 \(h^0\rightarrow e \mu \)

The latest experiment upper bound of decay \(h^0\rightarrow e \mu \) is smaller than \(0.035\%\) at \(95\%\) confidence level, which is detected by the CMS Collaboration. Al and \(A'l\) both appear in the non-diagonal terms of the slepton mass matrix. Considering the constraints from \(\mu \rightarrow e \gamma \) and \(\mu -e\) conversion in Au nuclei, we take \(M_{\tilde{E}}=A_E=\mu _{NL}=1\) TeV, \(A_{\tilde{E}}=\mu _Y=1.5\) TeV \(L_l=1\), \(L_F=0.006\), \(B_Y=0.94\) TeV, \(Y_{e4}=1.5\), \(Y_{e5}=0.8\), \(S_m^2=1 \mathrm{TeV}^2\), \(M_{Lf}^2=12000 \mathrm{GeV}^2\), \(m_1=m_2=0.5\) TeV, \(\mu =0.7\) TeV, \(\tan \beta =6\) and \(\tan \beta _{NL}=2\). The dotted (dashed, solid) line in Fig. 20 denotes the branching ratios of \(h^0\rightarrow e \mu \) versus \(A'l\) with \(Al=0.5(1,1.5)\) TeV. These three lines all increase quickly with the enlarging \(A'l\). So \(A'l\) play important roles to the numerical results. Although the larger Al, the smaller numerical results they are, the contributions from Al are very weak.

Fig. 20
figure 20

With \(Al=0.5(1,1.5)\) TeV, the branching ratios of \(h^0\rightarrow e \mu \) versus parameter \(A'l\) are plotted by the dotted (dashed, solid) line

At last, we discuss the effects from parameters \(Y_{e5}\) and \(M_{\tilde{E}}\). With \(m_1=m_2=0.5\) TeV, \(Al=1.5\) TeV, \(Y_{e4}=0.5\) and \(L_F=0.006\), the branching ratios varying with \(Y_{e5}\) are ploted in Fig. 21. The dotted, dashed and solid lines respectively correspond to \(M_{\tilde{E}}=1.1,1.4,1.7\) TeV. These three lines all slightly increase when \(Y_{e5}\) varies from 0.1 to 1.0. As \(Y_{e5}\) still increases from 1.0, the results have much more conspicuous enlargement. However, the total contributions from \(Y_{e5}\) are not so obvious. With the enlarging \(M_{\tilde{E}}\), the numerical results reduce more and more slowly.

5 Discussion and conclusion

We add exotic superfields \(\Phi _{NL}\), \(\varphi _{NL}\), Y and \(Y'\) to the BLMSSM, and this new model is named as the EBLMSSM. In \(\mathcal{W}_Y\), \(\lambda _4(\lambda _6)\) is the coupling coefficient of Y-lepton–exotic lepton and \(\tilde{Y}\)-lepton–exotic slepton. We assume \(\lambda _4^2=\lambda _6^2\) is a \(3\times 3\) squared matrix and its non-diagonal elements are related with the CLFV. Being different from the BLMSSM, the exotic slepton (sneutrino) of 4th and 5th generations mix together and form a \(4\times 4\) matrix. The Majorana particle, lepton neutralino \(\chi _L^0\), is corrected to be a \(5\times 5\) matrix due to the introduction of superpartners \(\psi _{\Phi _{NL}}\) and \(\psi _{\varphi _{NL}}\). The terms relating with exotic lepton (neutrino) and slepton (sneutrino) are also adjusted. In Sect. 3, we show the corresponding mass matrices and couplings of the EBLMSSM. The EBLMSSM has more abundant contents than that BLMSSM for the lepton physics.

Fig. 21
figure 21

With \(M_{\tilde{E}}=1.1(1.4,1.7)\) TeV, the branching ratios of \(h^0\rightarrow e \mu \) versus parameter \(Y_{e5}\) are plotted by the dotted, dashed and solid lines respectively

Considering the constraints from decays \(h^0\rightarrow \gamma \gamma \) and \(h^0\rightarrow VV,V=(Z,W)\), we study the CLFV processes \(l_j\rightarrow l_i \gamma \), muon conversion to electron in Au nuclei and the \(\tau \) decays in the framework of the EBLMSSM. Parameters \(Y_{e5}\) and \(M_{Lf}\) affect the numerical results in a certain degree. As the new introduced parameters in the EBLMSSM, \(\mu _Y,\tan {\beta }_{NL}, A_{\tilde{E}},A_E,L_l,M_{\tilde{E}}\) and \(v_N\) play important roles. Especially parameters \(L_f\) and \(L_F\) are all very sensitive parameters, which influence the numerical results very remarkably. Figures 13, 14 and 16 indicate that the enlarging \(L_f\) and \(L_F\) can easily improve the numerical results. Then, the 125.1 GeV Higgs boson decays with CLFV \(h^0\rightarrow l_il_j\) are discussed. As an important constraint, \(m_{h^0}=125.1\) GeV is regarded as an input parameter. Taking into account the constraints from the parameter space of decays \(l_j\rightarrow l_i \gamma \), muon conversion to electron in Au nuclei and the \(\tau \) decays, we analyze the numerical results for \(h^0\rightarrow l_il_j\) in EBLMSSM. Parameters \(\mu \) and \(A'l\) affect the CLFV processes in a certain degree. The effects from \(\tan \beta \) are very obvious. So \(\tan \beta \) is a sensitive parameter. Above all, due to the new particles introduced in the EBLMSSM, the numerical results can easily approach to the present experiment upper bounds.