Lepton flavor violating Higgs boson decays in seesaw models: new discussions

TL;DR

This paper analyzes lepton flavor violating Higgs decays in seesaw models at one-loop level, providing detailed calculations and comparisons for branching ratios across different models and neutrino mass scales.

Contribution

It offers new analytic expressions for LFVHD branching ratios in seesaw models and compares results across gauges and models, revealing relations between MSS and ISS predictions.

Findings

Maximal branching ratio around 10^{-11} in minimal seesaw.

Branching ratio increases with neutrino mass scale m_{n_6}.

Ratio of LFVHD rates in ISS and MSS models depends on m_{n_6} and x_X.

Abstract

The lepton flavor violating decay of the Standard Model-like Higgs boson (LFVHD), h->\mu\tau, is discussed in seesaw models at the one-loop level. Based on particular analytic expressions of Passarino-Veltman functions, the two unitary and 't Hooft Feynman gauges are used to compute the branching ratio of LFVHD and compare with results reported recently. In the minimal seesaw (MSS) model, the branching ratio was investigated in the whole valid range 10^{-9}-10^{15} GeV of new neutrino mass scale m_{n_6}. Using the Casas-Ibarra parameterization, this branching ratio enhances with large and increasing m_{n_6}. But the maximal value can reach only order of 10^{-11}. Interesting relations of LFVHD predicted by the MSS and inverse seesaw (ISS) model are discussed. The ratio between two LFVHD branching ratios predicted by the ISS and MSS is simply m^2_{n_6}\mu^{-2}_X, where \mu_X is the small neutrino mass scale in the ISS. The consistence between different calculations is shown precisely from analytical approach.

Figures

Tables

Table 1. Table 1: Couplings relating with LFVHD in seesaw models. Here, C i ​ j = ∑ c = 1 3 U c ​ i ν ​ U c ​ j ν ⁣ ∗ subscript 𝐶 𝑖 𝑗 superscript subscript 𝑐 1 3 subscript superscript 𝑈 𝜈 𝑐 𝑖 subscript superscript 𝑈 𝜈 𝑐 𝑗 C_{ij}=\sum_{c=1}^{3}U^{\nu}_{ci}U^{\nu*}_{cj} . The p 0 , p + subscript 𝑝 0 subscript 𝑝 p_{0},\,p_{+} and p − subscript 𝑝 p_{-} are incoming momenta of h ℎ h , G W + subscript superscript 𝐺 𝑊 G^{+}_{W} and G W − subscript superscript 𝐺 𝑊 G^{-}_{W} , respectively.

Vertex	coupling	Vertex	coupling
$h{W}^{+\mu }{W}^{-\nu }$	$ig{m}_W{g}_{\mu \nu }$	$h{G}_W^{+}{G}_W^{-}$	$-\frac{ig{m}_h^2}{2m_W}$
$h{G}_W^{+}{W}^{-\mu }$	$\frac{ig}{2}{\left(p_{+}-p_0\right)}_{\mu }$	$h{G}_W^{-}{W}^{+\mu }$	$\frac{ig}{2}{\left(p_0-p_{-}\right)}_{\mu }$
$\overline{n_i}{e}_a{W}_{\mu }^{+}$	$\frac{ig}{\sqrt{2}}{U}_{ai}^{\nu }{\gamma }^{\mu }{P}_L$	$\overline{e_a}{n}_i{W}_{\mu }^{-}$	$\frac{ig}{\sqrt{2}}{U}_{ai}^{\nu \ast }{\gamma }^{\mu }{P}_L$
$\overline{n_i}{e}_a{G}_W^{+}$	$-\frac{ig}{\sqrt{2}{m}_W}{U}_{ai}^{\nu }\left(m_{e_a}{P}_R-m_{n_i}{P}_L\right)$	$\overline{e_a}{n}_i{G}_W^{-}$	$-\frac{ig}{\sqrt{2}{m}_W}{U}_{ai}^{\nu \ast }\left(m_{e_a}{P}_L-m_{n_i}{P}_R\right)$
$h\overline{n_i}{n}_j$	$\frac{-ig}{2m_W}[C_{ij}(P_L{m}_{n_i}+P_R{m}_{n_j})$	$h\overline{e_a}{e}_a$	$-\frac{ig{m}_{e_a}}{2m_W}$
	$+C_{ij}^{\ast }(P_L{m}_{n_j}+P_R{m}_{n_i})]$

Equations

$$-\Delta \mathcal{L}=Y_{\nu ,aI}\overline{{\psi }_{L,a}}\tilde{\varphi }{N}_{R,I}+\frac{1}{2}\overline{(N_{R,I}{)}^c}{m}_{M,IJ}{N}_{R,J}+h.c.,$$

$$\displaystyle-\mathcal{L}^{\nu}_{\mathrm{mass}}\equiv\frac{1}{2}\overline{\nu^{\prime}_{L}}M^{\nu}(\nu^{\prime}_{L})^{c}+\mathrm{h.c.}=\frac{1}{2}\overline{\nu^{\prime}_{L}}\left(\begin{array}[]{cc}0&M_{D}\\ M_{D}^{T}&M_{N}\\ \end{array}\right)(\nu^{\prime}_{L})^{c}+\mathrm{h.c.},$$

$$U^{\nu T}{M}^{\nu }{U}^{\nu }={\hat{M}}^{\nu }=\mathrm{diagonal}(m_{n_1},m_{n_2},m_{n_3},m_{n_4},...,m_{n_{(K+3)}}),$$

$${\nu }_L^{\prime }=U^{\nu \ast }{n}_L,\mathrm{and}({\nu }_L^{\prime }{)}^c=U^{\nu }(n_L{)}^c,$$

$$P_L{\nu }_i^{\prime }={\nu }_{L,i}^{\prime }=U_{ij}^{\nu \ast }{n}_{L,j},\mathrm{and}{P}_R{\nu }_i^{\prime }={\nu }_{R,i}^{\prime }=U_{ij}^{\nu }{n}_{R,j},i,j=1,2,...,K+3,$$

$${\mathcal{L}}_{\mathrm{kin}}^{\mathrm{lep}}=i\overline{{\psi }_{L,a}}{\gamma }^{\mu }{D}_{\mu }{\psi }_{L,a}$$

$$-{\mathcal{L}}_{\mathrm{Y}}^{\mathrm{lep}}$$

$$M_{a(I+3)}^{\nu }{U}_{ai}^{\nu }{U}_{(I+3)j}^{\nu }{P}_R+M_{a(I+3)}^{\nu \ast }{U}_{(I+3)i}^{\nu \ast }{U}_{aj}^{\nu \ast }{P}_L=\left(\sum _{a=1}^3{U}_{ai}^{\nu }{U}_{aj}^{\nu \ast }\right)\left(m_{n_i}{P}_L+m_{n_j}{P}_R\right),$$

$$M_{ab}^{\nu }$$

$$U^{\nu †}{U}^{\nu }$$

$$M_{a(I+3)}^{\nu }{U}_{ai}^{\nu }{U}_{(I+3)j}^{\nu }$$

$$Y_{\nu ,aI}\overline{e_{L,a}}{N}_{R,I}{G}_W^{-}$$

$$\Gamma (h\to \mu \tau )\equiv \Gamma (h\to {\mu }^{-}{\tau }^{+})+\Gamma (h\to {\mu }^{+}{\tau }^{-})\simeq \frac{m_h}{8\pi }\left(|{\Delta }_L{|}^2+|{\Delta }_R{|}^2\right),$$

$${\Delta }_L^{(a)}$$

$${\Delta }_R^{(a)}$$

$${\Delta }_L^{(b)}$$

$${\Delta }_R^{(b)}$$

$${\Delta }_L^{(c+d)}$$

$${\Delta }_R^{(c+d)}$$

$$\mathrm{Div}[{\Delta }_L^{(a)}]$$

$$\mathrm{Div}[{\Delta }_L^{(b)}]$$

$$\sum _{i,j=1}^{K+3}\sum _{c=1}^3{m}_{n_i}^2{U}_{ai}^{\nu \ast }{U}_{ci}^{\nu }{U}_{cj}^{\nu \ast }{U}_{bj}^{\nu }$$

$$F_L^{(1)}$$

$$F_R^{(1)}$$

$$F_L^{(2)}$$

$$F_R^{(2)}$$

$$F_L^{(3)}$$

$$F_L^{(4)}$$

$$F_R^{(4)}$$

$$F_L^{(5)}$$

$$F_R^{(5)}$$

$$F_L^{(6)}$$

$$F_R^{(6)}$$

$$F_L^{(7)}$$

$$F_L^{(9)}$$

$$F_L^{(8)}$$

$$F_R^{(8)}$$

$$F_L^{(10)}$$

$$F_R^{(10)}$$

$$C_0$$

$$B_0^{(12)}$$

$$B_1^{(1)}(M_0^2,M^2)$$

Full text

I Introduction

After the Higgs boson was observed by ATLAS and CMS higgsdicovery , the LFVHD has been searched experimentally exLFVh , where upper bounds for branching ratios (Brs) of the decays $h\rightarrow\mu\tau,e\tau$ are order of $\mathcal{O}(10^{-2})$. Signals of LFVHD at future colliders have been discussed, where sensitivities for detecting these channel decays are shown to be $10^{-5}$ in the near future LFVcol . Up to now, the lepton flavor violating (LFV) decays of the standard-model-like and new Higgs bosons have been investigated in many models beyond the standard model (SM) Apo ; Apo1 ; EArganda ; iseesaw ; e1612.09290 ; moress ; LFVgeneral ; SUSY ; THDL2 ; LFVHDUgauge ; KHHung . Among them, the MSS MSS is the simplest that can explain successfully the recent neutrino data. Naturally, the mixing between different flavor neutrinos leads to many LFV processes from loop corrections. But it predicts very suppressed branching ratios (Br) of LFV decays of charged leptons. Recent studies on the Br of LFVHD were also shown to be very small EArganda . In contrast, the ISS ISSmodel , another simple extension of the SM, predicts much larger values of LFV branching ratios, including those of LFVHD iseesaw ; e1612.09290 . In fact, the Br of LFVHD in the ISS were calculated in many different ways in order to guarantee the consistence of the LFVHD amplitudes.

We stress that understanding the mechanism for generating loop corrections to Brs of LFVHD in simple models like the MSS and ISS is very important for studying LFVHD processes in other complicated models. That is why LFVHD predicted by these two models were discussed in many works, for example Apo ; Apo1 ; EArganda ; iseesaw ; e1612.09290 ; moress . In the ISS, recent results in iseesaw showed that branching ratios of LFVHD increase with increasing values of very heavy neutrino masses when the Casas-Ibarra method ibarra was applied to formulating the Yukawa couplings of heavy neutrinos 111We thank Dr. E. Arganda for this comment. But the Brs are always constrained by upper bounds because of the perturbative limit of the Yukawa couplings. Using the mass insertion approximation, a recent study e1612.09290 also calculated the Br of LFVHD in the ISS model in both unitary and ’t Hooft Feynman, where previous results in iseesaw were confirmed to be well consistent in the region of parameters containing large new neutrino mass scale $m_{n_{6}}$. The above discussions indicate that although one-loop contributions in both MSS and ISS arise from the same set of Feynman diagrams, the two models predict very different Br values. The reason is the appearance of a small mass scale $\mu_{X}$ in the ISS, which gives tiny contributions to the heavy neutrino masses, but affects strongly on the neutrino mixing matrix. Hence there should exist simple relations between two expressions of Brs predicted by the two models. These interesting relations were not discussed previously, therefore will be focused in this work. We will show that if $m_{n_{6}}$ is large enough, the ratio between Brs of LFVHD of the ISS and MSS is order of $m^{2}_{n_{6}}\mu^{-2}_{X}$, enough to explain clearly the LFVHD difference between two models.

Regarding the MSS, LFVHD was discussed mainly in ranges of $10^{2}-10^{7}$ GeV Apo ; EArganda , while the valid range of the new neutrino mass scale is from $\mathcal{O}(10^{-9})$ GeV to $\mathcal{O}(10^{15})$ GeV. In addition, a good estimation made in Ref. Apo suggested that the Br may enhance with increasing masses of heavy neutrinos, even when the Casas-Ibarra parameterization is used. We note that this parameterization are now still widely used to investigate the signal of seesaw models at recent colliders SSbound . As a result, possibilities that large Brs of LFVHD may exist in ranges of new neutrino mass scales that were not mentioned previously. Therefore, studies the LFVHD in the whole valid range as well as new approaches to compare well-known results and confirm consistent analytic formulas for calculating Br of LFVHD in seesaw models are still interesting and necessary. These are main scopes of this work. In particular, in order to guarantee the stability of numerical results at very large values of $m_{n_{6}}$, LFVHD processes will be computed using analytic expressions of Passarino-Veltman functions (PV functions) given in ref. LFVHDUgauge . Using a mathematica code based on these functions, we found that it is much easier and more convenient to increase the precision than using available numerical packages such as Looptools looptool . This makes our calculation different from all previous works. In addition, the one-loop contributions to LFVHD in both unitary and ’t Hooft Feynman gauges will be constructed using notations in LFVHDUgauge . Then we cross-check the consistence between total amplitudes calculated in two gauges, and the ones established in previous works Apo ; EArganda ; iseesaw . A detailed checking divergence cancellation will be presented analytically. For the MSS, after showing that Br of LFVHD is suppressed with small $m_{n_{6}}$, we will pay attention mainly to the region with large $m_{n_{6}}$. To guarantee the consistence of our investigation on LFVHD in the MSS, the connection between analytic formulas of LFVHD amplitudes in the two models MSS and ISS will be discussed deeply. In this work, Yukawa couplings of new neutrinos are only investigated following the Casas-Ibarra parameterization ibarra . This parameterization was used to investigate independently LFVHD processes predicted by the MSS and ISS in Refs. EArganda ; iseesaw , where other important properties of LFVHD were presented in details.

Our work is arranged as follows. Sec. II establishes notations and couplings of a general seesaw model needed for studying LFVHD. In Sec. III, we construct LFVHD amplitudes in two unitary and ’t Hooft Feynman gauges using notations of PV functions given in LFVHDUgauge . Then we prove the divergent cancellation and the consistence between two expressions of the LFVHD amplitudes. In Sec. IV, we show the choice of parameterizing the neutrino mixing matrices. After that, the Brs of LFVHD are numerically investigated. We will focus on new results of LFVHD in the MSS, and interesting relations between the Brs predicted by two models MSS and ISS. Sec. V summarizes new results of this work.

II General formalism and couplings for LFVHD

The general seesaw model is different from the Standard Model (SM) by $K$ additional right-handed neutrinos, $N_{R,I}\sim(1,1,0)$ with $I=1,2,...,K$ numixing . The new Lagrangian part is

[TABLE]

where $a=1,2,3$; I,J=1,2,…,K; $\psi_{L,a}=(\nu_{L,a},e_{L,a})^{T}$ are $SU(2)_{L}$ lepton doublets and $(N_{R,I})^{c}=C\overline{N_{R,I}}^{T}$. The Higgs bosons are also doublets $\phi=(G_{W}^{+},(h+iG_{Z}+v)/\sqrt{2})^{T}$ and $\widetilde{\phi}=i\sigma_{2}\phi^{*}$. Each of them consists of three Goldstone bosons of $W^{\pm}$ and $Z$ bosons; a neutral CP-even Higgs boson $h$ and the vacuum expectation value (VEV), $\langle\phi\rangle=\frac{v}{\sqrt{2}}=174$ GeV ($v=246$ GeV). Notations for flavor states of active neutrinos are $\nu_{L}=(\nu_{L,1},\;\nu_{L,2},\;\nu_{L,3})^{T}$ and $(\nu_{L})^{c}\equiv((\nu_{L,1})^{c},\;(\nu_{L,2})^{c},\;(\nu_{L,3})^{c})^{T}$. Notations for new neutrinos are $N_{R}=(N_{R,1},\;N_{R,2},...,\;N_{R,K})^{T}$, and $(N_{R})^{c}=((N_{R,1})^{c},\;(N_{R,2})^{c},...,\;(N_{R,K})^{c})^{T}$. In the bases of the original neutrinos, $\nu^{\prime}_{L}\equiv(\nu_{L},\;(N_{R})^{c})^{T}$ and $(\nu^{\prime}_{L})^{c}=((\nu_{L})^{c},\;N_{R})^{T}$, the Lagrangian part (1) generates the following mass term for neutrinos,

[TABLE]

where $M_{N}$ is a symmetric and non-singular $K\times K$ matrix, and $M_{D}$ is a $3\times K$ matrix, $(M_{D})_{aI}=Y_{\nu,aI}\langle\phi\rangle$. The matrix $M^{\nu}$ is symmetric, therefore it can be diagonalized via $(K+3)\times(K+3)$ matrix, $U^{\nu}$, satisfying the unitary condition, $U^{\nu\dagger}U^{\nu}=I$. We define

[TABLE]

where $m_{n_{i}}$ ($i=1,2,...,K+3$) are mass eigenvalues of the $(K+3)$ mass eigenstates $n_{L,i}$, i.e. physical states of neutrinos. Three light active neutrinos are $n_{L,a}$ with $a=1,2,3$. The relation between the flavor and mass eigenstates are

[TABLE]

where $n_{L}\equiv(n_{L,1},n_{L,2},...,n_{L,K+3})^{T}$.

In calculation, we will use a general notation of four-component (Dirac) spinor, $n_{i}$ ($i=1,2,..,K+3$), for all active and exotic neutrinos. Specifically, a Majorana fermion $n_{i}$ is defined as $n_{i}\equiv(n_{L,i},\;(n_{L,i})^{c})^{T}=n^{c}_{i}=(n_{i})^{c}$. The chiral components are $n_{L,i}\equiv P_{L}n_{i}$ and $n_{R,i}\equiv P_{R}n_{i}=(n_{L,i})^{c}$, where $P_{L,R}=\frac{1\pm\gamma_{5}}{2}$ are chiral operators. The similar definitions for the original neutrino states are $\nu_{a}\equiv(\nu_{L,a},\;(\nu_{L,a})^{c})^{T}$, $N_{I}\equiv((N_{R,I})^{c},\;N_{R,I})^{T}$, and $\nu^{\prime}=(\nu,\,N)^{T}$. The relations in (6) are rewritten as follows,

[TABLE]

where more precise expressions are $\nu_{L,a}=P_{L}\nu^{\prime}_{a}=U^{\nu*}_{ai}n_{L,i}$, $(N_{R,I})^{c}=P_{L}\nu^{\prime}_{I+3}=U^{\nu*}_{(I+3)j}n_{L,j}$, $(\nu_{L,a})^{c}=P_{R}\nu^{\prime}_{a}=U^{\nu}_{ai}n_{R,i}$, and $N_{R,I}=P_{R}\nu^{\prime}_{I+3}=U^{\nu}_{(I+3)j}n_{R,j}$ ($I=1,2,3,..,K$).

As usual, the covariant derivative is $D_{\mu}=\partial_{\mu}-igT^{a}W^{a}-ig^{\prime}YB_{\mu}$. We emphasize that the signs in $D_{\mu}$ will result in signs of couplings $hG^{\pm}_{W}W^{\pm}$ and $\overline{e}_{a}\nu_{a}W^{-}$. Correspondingly, the lepton flavor violating (LFV) couplings of $W^{\pm}$ boson to leptons are,

[TABLE]

where $a=1,2,3$; and $j=1,2,...,K+3$.

The Yukawa couplings that contribute to LFVHD are

[TABLE]

Using $(M_{D})_{aI}=M^{\nu}_{a(I+3)}$, and $N_{R,I}=\nu^{\prime}_{R,(I+3)}$, the last line in (9) changes in to the new form, $\frac{1}{v}h\overline{n_{i}}\left[M^{\nu}_{a(I+3)}U^{\nu}_{ai}U^{\nu}_{(I+3)j}P_{R}+M^{\nu*}_{a(I+3)}U^{\nu*}_{(I+3)i}U^{\nu*}_{aj}P_{L}\right]n_{j}$. It can be proved that

[TABLE]

which was given in EArganda ; iseesaw . A proof is as follows, based on the following properties of $M^{\nu}$ and $U^{\nu}$ defined in Eqs. (4) and (5),

[TABLE]

The first term in the left hand side of Eq. (10) will change exactly into the second term in the right hand side of Eq. (10), after mediate steps of transformation, namely

[TABLE]

From (12), the second term in the left hand side of (10) can be derived easily, $M^{\nu*}_{a(I+3)}U^{\nu*}_{(I+3)i}U^{\nu*}_{aj}=\left[M^{\nu}_{a(I+3)}U^{\nu}_{aj}U^{\nu}_{(I+3)i}\right]^{*}=\left[U^{\nu}_{aj}U^{\nu*}_{ai}m_{n_{i}}\right]^{*}=U^{\nu}_{ai}U^{\nu*}_{aj}m_{n_{i}}$. Finally, the Feynman rule for the vertex (10) with two Majorana leptons $h\overline{n_{i}}n_{j}$ must be expressed in a symmetric form 222 We thank Dr. E. Arganda for showing us this point, namely $-\frac{g}{4m_{W}}\sum_{i,j}\overline{n_{i}}\left[\left(m_{n_{i}}C_{ij}+m_{n_{j}}C^{*}_{ij}\right)P_{L}+\left(m_{n_{j}}C_{ij}+m_{n_{i}}C^{*}_{ij}\right)P_{R}\right]n_{j},$ where $C_{ij}=\sum_{c=1}^{3}U^{\nu}_{ci}U^{\nu*}_{cj}$ Apo ; spinor .

The couplings relating with $G^{\pm}_{W}$ are proved the same way, namely

[TABLE]

The vertices relating to LFVHD are collected in Table 1.

We note that the coupling $hG^{+}_{W}G^{-}_{W}$ in Table 1 is consistent with that given in hgwgw ; e1612.09290 .

The effective Lagrangian of the LFVHD is written as $\mathcal{L}^{LFV}=h\left(\Delta_{L}\overline{\mu}P_{L}\tau+\Delta_{R}\overline{\mu}P_{R}\tau\right)+\mathrm{h.c.}$, where $\Delta_{L,R}$ are scalar factors arising from loop contributions. The partial decay width is

[TABLE]

where $m_{h}\gg m_{2},m_{3}$ and $m_{2},m_{3}$ being masses of muon and tau, respectively. The on-shell conditions for external momenta are $p^{2}_{a}=m_{a}^{2}$ ($a=2,3$) and $p_{h}^{2}\equiv(p_{2}+p_{3})^{2}=m^{2}_{h}$, $m_{h}=125$ GeV. Next, $\Delta_{L,R}$ with be calculated at one-loop level, in two gauges of unitary and ’t Hooft Feynman.

III Analytic amplitudes and divergence cancellation

III.1 Amplitude in the unitary gauge and divergence cancellation

In the unitary gauge, the Feynman diagrams for a decay $h\rightarrow e^{-}_{a}e^{+}_{b}$ $(a<b)$ are presented in Fig. 1.

The loop contributions are written as $\Delta_{L,R}=\Delta^{(a)}_{L,R}+\Delta^{(b)}_{L,R}+\Delta^{(c+d)}_{L,R}$, where the three terms come from private contributions of diagrams 1a), 1b), and sum of contributions from two diagrams c) and d), respectively. The analytic expressions of contributions from the three diagrams 1a), c), and d) can be derived directly from LFVHDUgauge , except the diagram 1b) containing the coupling $h\overline{n_{i}}n_{j}$. An analytic expression of $\Delta^{(b)}_{L,R}$ is derived in appendix C. We have used Form form to cross-check our results. In addition, the total $\Delta_{L,R}$ is consistent with the result calculated in the ’t Hooft Feynman gauge, as we will show later. Expressions of LFVHD contributions in the unitary gauge are

[TABLE]

and

[TABLE]

Regarding $\Delta^{(b)}_{L,R}$, the contributions from $B^{(1)}_{1}=B^{(1)}_{1}(m_{W}^{2},m^{2}_{n_{i}})$ and $B^{(2)}_{1}$ are zeros because, for example, $B^{(1)}_{1}$ contains a factor $\sum_{j}U^{\nu}_{bj}m_{n_{j}}U^{\nu}_{cj}=\left(U^{\nu*}\hat{M}^{\nu}U^{\nu\dagger}\right)^{*}_{bc}=M^{\nu*}_{bc}=0$.

Divergence cancellation in the total amplitude is explained as follows. From divergent parts of the PV functions in Appendix A, the divergent parts of $\Delta^{(a)}_{L}$ and $\Delta^{(b)}_{L}$ are

[TABLE]

where the unitary property of $U^{\nu}$ is used to cancel the second term of $\mathrm{Div}[\Delta^{(a)}_{L}]$, namely $\sum_{i=1}^{K+3}U^{\nu*}_{ai}U^{\nu}_{bi}=\left(U^{\nu}U^{\nu\dagger}\right)_{ab}=0$. The second term of $\mathrm{Div}[\Delta^{(b)}_{L}]$ vanishes because $\sum_{i}U^{\nu*}_{ai}U^{\nu*}_{ci}m_{n_{i}}=\left(U^{\nu*}\hat{M}_{\nu}U^{\nu\dagger}\right)_{ac}=M^{\nu}_{ac}=0$ with all $a,c=1,2,3$. We simplify the first term of $\mathrm{Div}[\Delta^{(b)}_{L}]$ based on the following equalities

[TABLE]

Similarly, we have $\sum_{i,j=1}^{K+3}\sum_{c=1}^{3}2m^{2}_{n_{j}}U^{\nu*}_{aj}U^{\nu}_{ci}U^{\nu*}_{cj}U^{\nu}_{bj}=\sum_{i=1}^{K+3}2m^{2}_{n_{i}}U^{\nu*}_{ai}U^{\nu}_{bi}$. Inserting these two results into $\mathrm{Div}[\Delta^{(b)}_{L}]$ will give $\mathrm{Div}[\Delta^{(b)}_{L}]+\mathrm{Div}[\Delta^{(a)}_{L}]=0$. With $\Delta^{(c+d)}_{L}$, the divergent parts of the two terms $m_{a}^{2}B^{(1)}_{1}$ and $m_{b}^{2}B^{(2)}_{1}$ vanish because of the GIM mechanism, while two sums $[B^{(1)}_{1}+B^{(2)}_{1}]$ and $[B^{(1)}_{0}-B^{(2)}_{0}]$ are finite. Hence, $\Delta_{L}$ is finite. $\Delta_{R}$ has the same conclusion.

III.2 Amplitude in the ’t Hooft Feynman gauge.

In the ’t Hooft Feynman gauge, there are ten form factors $F^{(i)}_{L,R}$, ($i=1,2,..,10$) corresponding to ten diagrams shown in Fig. 1 of Refs. EArganda ; iseesaw . The total contribution is $\Delta_{L,R}=\sum_{i=1}^{10}F^{i}_{L,R}$. Formulas of $F^{(i)}_{L,R}$ in terms of PV functions defined in LFVHDUgauge are as follows,

[TABLE]

where $B_{ai}=U^{\nu*}_{ai},B^{*}_{bj}=U^{\nu}_{bj}$, $C_{ij}=\sum_{c=1}^{3}U^{\nu}_{ci}U^{\nu*}_{cj}$, and $D=4-2\epsilon$ is the integral dimension defined in Appendix A. Although $F^{(7)}_{L,R}$ and $F^{(9)}_{L,R}$ contain $B$-functions, they are finite because of the GIM mechanism. Hence it can be replaced with $D=4$. Because $B^{(12)}_{0}=B^{(12)}_{0}(m_{W}^{2},m_{W}^{2})$ in $F^{(4)}_{R}$ and $F^{(5)}_{L}$ do not depend on $m_{n_{i}}$, therefore vanish because of the GIM mechanism. They will be ignored from now on.

Although our notations of PV functions are different from those in EArganda ; iseesaw , transformations between two sets of notations are, (see a detailed proving in Appendix B)

[TABLE]

The PV functions used in our work were checked to be consistent with Looptools looptool , see details in KHHung . The differences between our results and those shown in iseesaw are minus signs in $F^{(4)}_{L,R}$ and $F^{(5)}_{L,R}$. Our formulas are consistent with the results presented in Ref. e1612.09290 333The correct Feynman rule for the coupling $h\overline{n_{i}}n_{j}$ gives consistent $F^{(1,2)}_{L,R}$ with those in Ref. iseesaw ., where the authors confirmed that these signs do not affect the results given in Ref. iseesaw .

Now we will check the consistence between total amplitudes calculated in two gauges. Regarding to triangle diagrams with two internal neutrino lines, the deviation of contributions in two gauge are determined as follows,

[TABLE]

where useful equalities of B-functions are used bardin . In addition, $C_{ij}$ in the first line of (24) is simplified using the same trick given in (18). Similarly, other deviations are

[TABLE]

where $B_{0,1,2}\equiv B_{0,1,2}(m^{2}_{n_{i}},m^{2}_{W})$. Then, it can be seen easily that $\delta_{1}+\delta_{2}+\delta_{3}=0$. Hence, the total amplitudes calculated in two gauges are the same.

IV LFVHD in the minimal and inverse seesaw models

IV.1 Parameterization the neutrino mixing matrix

To start, we consider a general expression of the neutrino mixing matrix $U^{\nu}$ numixing ,

[TABLE]

where $\mathbf{O}$ is a $3\times K$ null matrix, $U$ and $V$ are $3\times 3$ and $K\times K$ unitary matrices, respectively. The $\Omega$ is a $(K+3)\times(K+3)$ unitary matrix that can be formally written as

[TABLE]

where $R$ is a $3\times K$ matrix where absolute values of al elements are smaller than unity. The unitary matrix $U=U_{\mathrm{PMNS}}$ is the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix upmns .

The mass matrices of neutrinos are written as follows,

[TABLE]

where $m_{n_{i}}$ is the physical masses of all neutrinos,

[TABLE]

and $c_{ab}\equiv\cos\theta_{ab}$, $s_{ab}\equiv\sin\theta_{ab}$. In the normal hierarchy scheme, the best-fit values of neutrino oscillation parameters are given as actnuUpdate 444 Updated neutrino data can be found in pdg2016 . But our main results are unchanged

[TABLE]

where $\Delta m^{2}_{a1}=m^{2}_{n_{a}}-m^{2}_{n_{1}}$ ($a=2,3$). In this work, other parameters will be fixed as $\delta=\alpha=\beta=0$.

The condition of seesaw mechanism for neutrino mass generation is $|M_{D}|\ll|M_{N}|$, where $|M_{D}|$ and $|M_{N}|$ denote characteristic scales of $M_{D}$ and $M_{N}$, resulting in useful relations 555We thank LE Duc Ninh for pointing out factors $1/2$ in the last relation in (34). numixing ,

[TABLE]

Based on the second relation in (34), the matrix $M_{D}$ can be parameterized via a general $K\times 3$ matrix $\xi$, which satisfies the only condition $\xi^{T}\xi=I_{3}$ ibarra ; numixing ; EArganda , namely

[TABLE]

where $U_{N}$ is an unitary matrix diagonalizing $M_{N}$, $U_{N}^{T}M_{N}U_{N}=M^{d}_{N}=\mathrm{diag}(M_{1},M_{2},...,M_{K})$.

In the MSS mentioned in Apo ; EArganda , the particle content is different from the Standard Model (SM) by three additional right-handed neutrinos ($K=3$), $N_{R,I}\sim(1,1,0)$ with $I=1,2,3$. New notations of neutrino mass matrices are $m_{D}\equiv M_{D}$, and $m_{M}\equiv M_{N}$. They are the respective $3\times 3$ Dirac and Majorana mass matrices corresponding to the first and second term of (1), $(m_{D})_{iJ}=Y_{\nu,iJ}\langle\phi\rangle$, and $(m_{M})_{iJ}=m_{M,iJ}$. The matrix $m_{M}$ is real, symmetric and non-singular.

The mixing matrix in the ISS model considered in ref. iseesaw can be found approximately using the above general discussion with $K=6$. Relations of notations between two parameterizations in iseesaw and numixing are

[TABLE]

where $\mathcal{O}$ is the $3\times 3$ matrix with all elements being zeros. From the definition of the inverse matrix, $M^{-1}_{N}M_{N}=M_{N}M^{-1}_{N}=I_{6}$, we derive that

[TABLE]

where $M$ is defined as $M=M_{R}\mu_{X}^{-1}M_{R}^{T}$ iseesaw . From (34), we then find that numixing

[TABLE]

These two expressions are consistent with those given in iseesaw ; numixing , giving a parameterization of $m_{D}$ as follows,

[TABLE]

where $U_{M}$ satisfies $M=U_{M}^{*}\mathrm{diag}(M_{1},\;M_{2},\;M_{3})U^{\dagger}_{M}$ and $\xi^{\prime}$ is a complex orthogonal matrix satisfying $\xi^{\prime}\xi^{\prime T}=I_{3}$. The mixing matrix $U^{\nu}$ now is a $9\times 9$ matrix.

In order to compare and mark relations between LFVHD in two MSS and ISS models, we will pay attention to only simply cases of choosing parameters. In the MSS model, the choice is $\xi=U_{N}=I_{3}$, leading to following simple expressions of Eqs. in (34), namely

[TABLE]

In the ISS model, from (43) we see that $m_{D}$ is parameterized in terms of many free parameters, hence it is enough to choose that $\mu_{X}=\mu_{X}I_{3}$. This parameter is a new scale making the most important difference between the neutrino mixing matrices in the ISS and MSS. We also assume that $M_{R}=\hat{M}_{R}=\mathrm{diag}(M_{R_{1}},\,M_{R_{2}},\,M_{R_{3}})$ and $\xi^{\prime}=I_{3}$. With $|\mu_{X}|\ll|M_{R}|$ we have

[TABLE]

We can see that both $\hat{M}_{R}$ (ISS) and $M_{N}$ (MSS) play roles as exotic neutrino mass scales. Therefore, they are identified as neutrino masses in both models, $\hat{M}_{R}=M_{N}=\mathrm{diag}(m_{n_{4}},\,m_{n_{5}},\,m_{n_{6}})$. The differences between two models now are two mixing matrix $V$ in (45) and $R$, and the $\mu_{X}$ scale, which does not appear in the MSS model. The $\mu_{X}$ plays special roles in the ISS model via its appearance in the second sub-matrix of the mixing matrix $R$ given in (42). A simple relation between largest elements of $R$ matrices in two models is

[TABLE]

where $m_{n_{6}}$ now is considered as exotic neutrino mass scale, $m_{n_{4}}\leq m_{n_{5}}\leq m_{n_{6}}$. The relation (46) is the main reason that explains why the Br of LFVHD predicted from the ISS is much larger than that from the MSS.

In the following, we will discuss on LFVHD in the MSS model. The results of LFVHD in the ISS model can be derived from discussion in the MSS model based on (46).

IV.2 Discussion on LFVHD

In the MSS model, our investigation will use three physical masses of exotic neutrinos, $m_{n_{4,5,6}}$, as free parameters. The matrix $M_{D}$ can be derived from relations (35), i.e $M_{D}=iU^{*}_{\mathrm{PMNS}}\left(M^{d}_{N}\,\hat{m}_{\nu}\right)^{1/2}$. As a result, the mixing matrix $U^{\nu}$ is written as a function of physical neutrino masses and $U_{\mathrm{PMNS}}$. To determine constrains of heavy neutrino masses $m_{n_{6}}$, we base on relations in (34), which suggest that $m_{n_{6}}\times m_{n_{3}}\simeq|M_{D}|^{2}<6\pi\times 174^{2}$, because of the perturbative limit of the Yukawa couplings $Y_{\nu,ij}$ iseesaw . Combing with the active neutrino data given in (33), where at least one active neutrino mass is not smaller than $\sqrt{\Delta m_{31}^{2}}=5\times 10^{-11}$ GeV, we get an upper constrain, $m_{n_{6}}<8\times 10^{15}$ GeV, when $m_{n_{1}}\ll\sqrt{\Delta m_{31}^{2}}$. The lower constrain is $m_{n_{6}}>|M_{D}|>m_{n_{3}}>5\times 10^{-11}$ GeV. Numerical illustrations are shown in Fig. 2, where three heavy neutrino masses are non-degenerate, $3m_{n_{4}}=2m_{n_{5}}=m_{n_{6}}$, and $m_{n_{1}}=10^{-12}\,\mathrm{GeV}\ll\sqrt{\Delta m_{31}^{2}}$.

The left panel of Fig. 2 presents Br$(h\rightarrow e_{a}e_{b})$ as functions of $m_{n_{6}}$. Unlike previous works such as Apo ; EArganda , heavy neutrinos masses were not considered at the interesting scale above $10^{10}$ GeV, where leptogenesis can be successful explained in the MSS frame work leptoMSS . More important, large values of heavy neutrinos may give large Br of LFVHD, as we have seen numerically. Unfortunately, values of $m_{n_{6}}\leq 8\times 10^{15}$ GeV gives an upper bound Br$(h\rightarrow\mu\tau)\leq\mathcal{O}(10^{-11})$. For other two decays, we get the relations Br$(h\rightarrow e\tau)\simeq$Br$(h\rightarrow\mu\tau)=(m^{2}_{\tau}/m^{2}_{\mu})$Br$(h\rightarrow e\mu)\simeq 287\times\mathrm{Br}(h\rightarrow e\mu)$. Hence, we just focus on the Br$(h\rightarrow\mu\tau)$.

The right panel of Fig. 2 shows values of Br$(h\rightarrow\mu\tau)$ in the whole valid range of $m_{n_{6}}$, namely $10^{-10}<m_{n_{6}}<8\times 10^{15}$ [GeV], where $U^{\nu}$ is considered up to $\mathcal{O}(R^{k})$. Each curve separates into three different parts. In the part with very heavy exotic neutrino masses, $m^{2}_{n_{6}}\gg m^{2}_{h},m^{2}_{W}$, i.e. $m_{n_{6}}>\mathcal{O}\left(10^{4}\right)$, we found a simple relation: Br$(h\rightarrow\mu\tau)=6.3\times 10^{-44}m_{n_{6}}^{2}\,[\mathrm{GeV}^{2}]$. On the other hand, for the part with very small exotic neutrino masses, $m^{2}_{n_{6}}\ll m^{2}_{\mu},m^{2}_{\tau}$, i.e. $m_{n_{6}}<\mathcal{O}\left(10^{-3}\right)$, there appears a new relation: Br$(h\rightarrow\mu\tau)=\dfrac{8.7\times 10^{-52}}{(m_{n_{6}}\mathrm{GeV})^{4}}$, when the matrix $\Omega$ is calculated up to $\mathcal{O}(R^{2})$. This will lead to the maximal values of Br$(h\rightarrow\mu\tau)\leq 10^{-11}$, the same order with large $m_{n_{6}}\sim\mathcal{O}(10^{15})$ GeV. If the matrix $\Omega$ is calculated more exactly, the Br$(h\rightarrow\mu\tau)$ will decrease significantly with small $m_{n_{6}}$, but will not change with large $m_{n_{6}}$. This can be explained from the conditions of the matrix $\Omega$, which is written in terms of the power series in $R$. If $m_{n_{6}}$ is small, $R\sim\sqrt{|m_{\nu}|/m_{n_{6}}}$ will be large as $m_{n_{6}}\rightarrow|M_{D}|\rightarrow|m_{\nu}|$. The calculation will be less accurate with smaller power $k$ included in $\Omega$. We consider more cases of $U^{\nu}$ where the matrix $\Omega$ in (27) is considered up to order $\mathcal{O}(R^{8})$. We conclude that the Br$(h\rightarrow\mu\tau)$ is very suppressed with small masses of exotic neutrinos. In contrast, large $m_{n_{6}}$ results in $|R|\ll 1$. Therefore, it is enough to consider the mixing matrix $U^{\nu}$ with order of $\mathcal{O}(R^{2})$ in the region where $m_{n_{6}}\geq 0.1$ GeV. In conclusion, to find large Br$(h\rightarrow\mu\tau)$, we just consider the region with large $m_{n_{6}}$.

To explain why large Br$(h\rightarrow\mu\tau)$ corresponds to large $m_{n_{6}}$, we pay attention to the properties of the mixing matrix $U^{\nu}$, the PV-functions and factors relating with them in the expressions of $\Delta^{(a)}_{L,R}$, $\Delta^{(b)}_{L,R}$, and $\Delta^{(c+d)}_{L,R}$. When $m^{2}_{n_{I}}\gg m^{2}_{h},m^{2}_{W}$, the terms with factors $m^{2}_{n_{I}}$ will give dominant contributions. The PV functions containing $m^{2}_{n_{I}}$ will have the following properties: $B_{0,1,2}(m^{2}_{n_{i}})=\mathcal{O}(10)$, $C_{0,1,2}(m^{2}_{n_{i}})\sim\ln(m^{2}_{n_{6}})/m^{2}_{n_{6}}$. Hence the largest contributions will come from $m^{2}_{n_{6}}B_{0,1,2}\sim m^{2}_{n_{6}}$ in $\Delta^{(a+c+d,b)}_{L,R}$ and $m^{4}_{n_{6}}C_{0,1,2}\sim[\ln\,m^{2}_{n_{6}}]m^{2}_{n_{6}}$ in $\Delta^{(b)}_{L,R}$.

The largest component of the matrix $R$ satisfies $R\sim\mathcal{O}\left(\sqrt{\frac{|\hat{m}_{\nu}|}{m_{n_{6}}}}\right)$. As a result, the mixing matrix elements in $\Delta^{(a+c+d)}_{L,R}$ and $\Delta^{(b)}_{L,R}$ will results in the following factors: $U^{\nu*}_{a(I+3)}U^{\nu}_{b(I+3)}=|R_{aI}|^{2}\sim\frac{|\hat{m}_{\nu}|}{m_{n_{6}}}$. There are new factors in the $\Delta^{(b)}_{L,R}$: $U^{\nu*}_{a(I+3)}U^{\nu}_{c(I+3)}U^{\nu*}_{c(J+3)}U^{\nu*}_{b(J+3)}\sim\frac{|\hat{m}_{\nu}|^{2}}{m^{2}_{n_{6}}}$. Hence the largest contribution to the total gives $\Delta_{L,R}\sim m_{n_{6}}$ with very large $m_{n_{6}}$, implying Br$(h\rightarrow\mu\tau)\sim m^{2}_{n_{6}}$. The correlations between terms with and without factors $m^{2}_{n_{i}}$ are shown in the Fig. 3. Terms without factors $m_{n_{i}}^{2}$ are dominant with tiny $m_{n_{6}}$ but they are very suppressed with large $m_{n_{6}}$.

The above discussions lead to new interesting results for LFVHD predicted by the MSS model, which were not concerned previously: i) the Br can reach values of order $10^{-11}$ with large values of heavy neutrino masses satisfying the perurbative limit; ii) the Br enhances with increasing $m_{n_{6}}$ above $10^{5}$ GeV. In addition, the maximal Br$(h\rightarrow\mu\tau)$ reaches the values of $10^{-33}-10^{-32}$ with $m_{n_{6}}\in[10^{2},10^{4}]$ GeV. We will show the relation between these interesting values and maximal values of Br$(h\rightarrow\mu\tau)$ predicted by the ISS.

We realize that the property of Br$(h\rightarrow\mu\tau)\sim m^{2}_{n_{6}}$ agrees very well with the approximate expression shown in Apo . In particular, Br$(h\rightarrow\mu\tau)\sim m^{4}_{n_{6}}\times|F_{N}|^{2}$, where $F_{N}\sim R^{2}\sim m^{-1}_{n_{6}}$ relating with active-heavy neutrino mixing elements in $U^{\nu}$. We believe that large values of the Br predicted in Apo arise from the reason that recent neutrino oscillation data could not be applied at that times. The numerical values of $F_{N}$ chosen in Apo may keep large contributions that should vanish because of the GIM mechanism.

Although the maximal Br of LFVHD predicted by the MSS is much smaller than the prediction from the ISS model given in EArganda ; iseesaw , the behave of the curve presenting Br$(h\rightarrow\mu\tau)$ shown in Fig. 3 have the same form with Br$(h\rightarrow\mu\tau)$ calculated in the ISS. The reason is as follows. If the exotic neutrino masses are fixed the same values in the two models, $m_{M}=M_{R}=\mathrm{diag}(m_{n_{4}},m_{n_{5}},m_{n_{6}})$, the important quantity making different contributions to LFVHD is the parametrization of $m_{D}$, see two Eqs. (35) and (43) for the MSS and ISS, respectively.

This leads to the different structures of the $R$ matrices. The largest components of $R$ in the MSS are $R^{\mathrm{MSS}}_{aI}\sim\sqrt{\frac{|\hat{m}_{\nu}|}{|m_{n_{6}}|}}$ with $I>3$, while those in the ISS are $R^{\mathrm{ISS}}_{aI}\sim\sqrt{\frac{|\hat{m}_{\nu}|}{\mu_{X}}}$. Hence, in general the ISS mixing factors are larger than those of MSS a common factor $\sqrt{\frac{|\hat{m}_{n_{6}}|}{|\mu_{X}|}}$. It makes the prediction of Br of LFVHD by the ISS be much larger than the prediction by the MSS, provided large $m_{n_{6}}$ but small $\mu_{X}$. Unlike the MSS, where mass scale $m_{n_{6}}$ can be as large as $\mathcal{O}(10^{15})$ GeV, values of $m_{n_{6}}$ in the ISS are constrained by relation (42), i.e. $m^{2}_{n_{6}}|\hat{m}_{\nu}|/\mu_{X}=|m_{D}|^{2}<174^{2}\times 6\pi\,[\mathrm{GeV}]^{2}$. Hence, small $\mu_{X}$ will give small upper bounds of $m_{n_{6}}$, and large Br$(h\rightarrow\mu\tau)$ will depend complicatedly on these two parameters. The left panel of Fig. 4 shows possible values of Br$(h\rightarrow\mu\tau)$ in the allowed regions of $\mu_{{}_{X}}$ and $m_{n_{6}}$. Our numerical results are well consistent with previous work iseesaw . In addition, by adding a factor $\sqrt{\frac{|m_{n_{6}}|}{|\mu_{X}|}}$ into $R^{\mathrm{MSS}}$ and using the analytic expressions of $\Delta^{\mathrm{MSS}}_{L,R}$ we get a very consistent results of Br$(h\rightarrow\mu\tau)$ predicted by the ISS, see an illustration in the left panel of Fig. 4. This confirms again the consistence of our calculation for LFVHD in the MSS and ISS.

There is an interesting relation between two LFVHD amplitudes calculated in the two models, as drawn in the right panel of Fig. 4. Here, $|\Delta^{\mathrm{ISS}}_{R}|\mu_{X}m^{-2}_{n_{6}}$ and $|\Delta^{\mathrm{MSS}}_{R}|m^{-1}_{n_{6}}$ are considered as functions of $m_{n_{6}}$. We have checked numerically that $|\Delta^{\mathrm{ISS}}_{R}|\mu_{X}m^{-2}_{n_{6}}$ does not depend on $\mu_{X}$, and consistent with conclusion in iseesaw . It can be seen as follows. The dependence of $m_{D}$ and $R^{\mathrm{ISS}}$ on $M_{R}$ and $\mu_{X}$ can be separate into two parts. The first is the correlation between elements of these matrices in order to give correct experimental values of active neutrino data. And the second is the simple dependence on the scales of $m_{n_{6}}$ and $\mu_{X}$. In the ISS, $R^{\mathrm{ISS}}_{aI}=U^{\nu}_{a(I+3)}\sim\mu_{X}^{-1/2}$ and do not depend on $m_{n_{6}}$. Now, if we pay attention to the region with large $m_{n_{6}}$, the terms like $m_{n_{i}}^{2}B_{0,1,2}$ are dominant contributions to $\Delta_{L,R}$ because of the factors $m_{n_{i}}^{2}$. As a result, $\Delta^{(a+c+d)}_{L,R}$ containing a factor $U^{\nu*}_{ai}U^{\nu}_{bi}\sim\mu^{-1}_{X}$ will give an overall factor $\mu_{X}^{-1}m^{2}_{n_{6}}$. Hence $\Delta^{(a+c+d)}_{L,R}\mu_{X}m^{-2}_{n_{6}}$ may be constant, following the property of $B$-functions. On the other hand, $\Delta^{(b)}_{L,R}$ contains $U^{\nu*}_{ai}U^{\nu*}_{cj}U^{\nu}_{ci}U^{\nu}_{bj}\sim\mu^{-1}_{X}$ or $\mu^{-2}_{X}$, depending on both indices $i$ and $j$ or only one larger than 3. Because both $\Delta^{(a+c+d)}_{L,R}$ and $\Delta^{(b)}_{L,R}$ are still divergent, terms with $\mu^{-2}_{X}$ must vanish in order to guarantee a finite $\Delta_{L,R}$. This results in a common factor $\mu_{X}^{-1}m^{2}_{n_{6}}$ for $\Delta_{L,R}$. In the right panel of Fig. 4, values of $\mu_{X}m^{-2}_{n_{6}}\Delta^{(a+c+d)}_{L,R}$ and $\mu_{X}m^{-2}_{n_{6}}\Delta^{(b)}_{L,R}$ correspond to $\Delta_{\epsilon}=0$. But we checked numerically that $\mu_{X}^{-1}m^{2}_{n_{6}}\Delta_{L,R}$ is independent with $\Delta_{\epsilon}$. In addition, we can see that $\mu_{X}m^{-2}_{n_{6}}\Delta^{(a+c+d)}_{L,R}$ and $\mu_{X}m^{-2}_{n_{6}}\Delta^{(b)}_{L,R}$ always have opposite signs, which is consistent with the fact that divergences contained in them are really canceled. Two absolute contributions from $\Delta^{(a+c+d)}_{L,R}$ and $\Delta^{(b)}_{L,R}$ are the same order, and nearly degenerate with large $m_{n_{6}}$. They start canceling strongly each other from the electroweak range of $m_{n_{6}}$, giving a very small $\mu_{X}m^{-2}_{n_{6}}\Delta_{L,R}$. It is $10^{-5}$ times smaller than values of $\mu_{X}m^{-2}_{n_{6}}\Delta^{(b)}_{L,R}$.

The above discussion is the same for both models ISS and MSS, where $m^{-1}_{n_{6}}\Delta_{L,R}$ is the function considered in the MSS. The numerical results are also shown in the right panel of the Fig. 4. Consider a region $10\leq m_{n_{6}}\leq 10^{4}$ GeV, there is an equality that $m^{-1}_{n_{6}}\Delta^{\mathrm{MSS}}_{L,R}=\Delta^{\mathrm{ISS}}_{L,R}\mu_{X}m^{-2}_{n_{6}}$, implying $\mathrm{Br^{ISS}}(h\rightarrow\mu\tau)=\frac{m^{2}_{n_{6}}}{\mu^{2}_{X}}\mathrm{Br^{MSS}}(h\rightarrow\mu\tau)$. From previous discussion, where $\mathrm{Br^{MSS}}(h\rightarrow\mu\tau)\leq 10^{-32}$, we can derive the maximal $\mathrm{Br^{ISS}}(h\rightarrow\mu\tau)\leq 10^{-32}\times\mathcal{O}((10^{4}/10^{-9})^{2})=\mathcal{O}(10^{-6})$.

We can also estimate the maximal value of Br$(h\rightarrow\mu\tau)$ based on the numerical result shown in Fig. 4. If $m_{n_{6}}\geq 10^{5}$ GeV, we have $\Delta_{R}\simeq 10^{-24}\mu_{X}^{-1}m^{2}_{n_{6}}$, where small $\Delta_{L}$ is ignored. Equivalently, we have Br$(h\rightarrow\mu\tau)\simeq 10^{-45}\mu_{X}^{-2}m^{4}_{n_{6}}$. The condition of perturbative limit gives $m^{2}_{n_{6}}\times 5\times 10^{-11}/\mu_{X}=|m_{D}|^{2}\leq 174^{2}\times 6\pi$, leading to $\mu_{X}^{-2}m^{4}_{n_{6}}\leq\mathcal{O}(10^{36})$. Hence in the region of lagre $m_{n_{6}}\geq 10^{5}$ GeV, Br$(h\rightarrow\mu\tau)$ can reach maximal value of $\mathcal{O}(10^{-9})$. If $m_{n_{6}}<10^{5}$ GeV, the allowed region in the left panel of Fig. 4 shows that Br$(h\rightarrow\mu\tau)$ can reach values of $\mathcal{O}(10^{-7})$ only if $m_{n_{6}}$ is few TeV, $\mu_{X}$ is order of $10^{-9}$ GeV, and $m_{D}$ gets values very close to the perturbative limit.

V Conclusion

In this work, the LFVHD in the MSS and ISS models have been discussed where we have focused on new aspects that were not shown in previous works. We calculated the amplitude of the LFVHD using new analytical expressions of PV-functions discussed recently. From this we have checked the consistence of our results in many different ways: comparing them with results of previous works, calculating in two gauges of unitary and ’t Hooft-Feynman, checking analytically the divergent cancellation of the total amplitude. In the MSS framework, we investigated numerically the Br$(h\rightarrow\mu\tau)$ in the valid and large range of exotic neutrino mass scale, from $10^{-10}$ GeV to $10^{16}$ GeV. When applying the Casas-Ibarra parameterization to Yukwa couplings of heavy neutrinos, we found a new result that Br$(h\rightarrow\mu\tau)\sim m^{2}_{n_{6}}$ with large $m_{n_{6}}$, because the mixing matrix elements affecting mostly the LFVHD amplitude by factors of $m_{n_{6}}^{-1/2}$. But in the valid region of perturbative requiring $m_{n_{6}}<10^{16}$ GeV, the Br$(h\rightarrow\mu\tau)$ reaches maximal values of $\mathcal{O}(10^{-11})$, still far from the recent experimental consideration. Anyway, this may be a hint to improve the MSS to more relevant models predicting higher values of Br$(h\rightarrow\mu\tau)$, for example the ISS. In this model, the largest mixing factors contributing to LFVHD amplitude do not depend on the exotic neutrino mass scale $m_{n_{6}}$ but consist of a factor $\mu_{X}^{-1}$. Hence, if two models have the same neutrino mass scale, and the neutrino mixing matrices obey the Casas-Ibarra parameterization, there will be a very simple relation that $\mathrm{BR^{ISS}}(h\rightarrow\mu\tau)/\mathrm{BR^{MSS}}(h\rightarrow\mu\tau)\simeq m^{2}_{n_{6}}\mu^{-2}_{X}$. This explains why the signal of LFVHD in the ISS is extremely significant than that in MSS. But the perturbative condition does not allow both large $m_{n_{6}}$ and small $\mu_{X}$, which can predict large Br$(h\rightarrow\mu\tau)$. Hence, maximal Br$(h\rightarrow\mu\tau)$ is still $\mathcal{O}(10^{-7})$ with few TeV of heavy neutrino mass scale. Our discussion on LFVHD of the MSS suggests that Br$(h\rightarrow\mu\tau)$ may be large in the extended versions of the MSS which allow very large $m_{n_{6}}$. Finally, although we presented here a different way to calculate the LFVHD, our numerical results for the ISS are well consistent with those noted in previous works iseesaw ; e1612.09290 .

Acknowledgments

LTH thanks Professor Thomas Hahn for useful discussions on Looptools. We are especially thankful Dr. Ernesto Arganda and authors of Refs. EArganda ; iseesaw ; e1612.09290 for important comments to correct our calculations and help us understand more deeply the LFVHD in the ISS. We thank Dr. Julien Baglio and Dr. LE Duc Ninh for helpful discussions. This research is funded by Ministry of Education and Training under grant number B2017_SP2_06.

Appendix A One loop Passarino-Veltman functions

Calculation in this section relates with one-loop diagrams in the Fig. 1. The analytic expressions of the PV-functions are given in LFVHDUgauge and they were derived from the general forms given in Hooft , using only the conditions of very small masses of tau and muon. They are consistent with bardin . The denominators of the propagators are denoted as $D_{0}=k^{2}-M_{0}^{2}+i\delta$, $D_{1}=(k-p_{1})^{2}-M_{1}^{2}+i\delta$ and $D_{2}=(k+p_{2})^{2}-M_{2}^{2}+i\delta$, where $\delta$ is infinitesimally a positive real quantity. The scalar integrals are defined as

[TABLE]

where $i=1,2$. In addition, $D=4-2\epsilon\leq 4$ is the dimension of the integral; $M_{0},~{}M_{1},~{}M_{2}$ are masses of virtual particles in the loop. The momenta satisfy conditions: $p^{2}_{1}=m^{2}_{1},~{}p^{2}_{2}=m^{2}_{2}$ and $(p_{1}+p_{2})^{2}=m^{2}_{h}$. In this work, $m_{h}$ is the SM-like Higgs mass, $m_{1,2}$ are lepton masses. The tensor integrals are

[TABLE]

The PV functions are $B^{(i)}_{0,1}$, $B^{(12)}_{0}$ and $C_{0,1,2}$. The functions $C_{0,1,2}$ are finite while the remains are divergent. We define the common divergent part as $\Delta_{\epsilon}\equiv\frac{1}{\epsilon}+\ln 4\pi-\gamma_{E}+\ln\mu^{2}$ where $\gamma_{E}$ is the Euler constant. Then the divergent parts of the above scalar factors are $\mathrm{Div}[B^{(i)}_{0}]=\mathrm{Div}[B^{(12)}_{0}]=\Delta_{\epsilon}$, and $\mathrm{Div}[B^{(1)}_{1}]=-\mathrm{Div}[B^{(2)}_{1}]=\frac{1}{2}\Delta_{\epsilon}$.

For simplicity in calculation we use approximative forms of PV functions where $p_{1}^{2},p_{2}^{2}\rightarrow 0$. The function $C_{0}$ was given in LFVHDUgauge consistent with that discussed on bardin , namely

[TABLE]

where $R_{0}(x_{0},x_{i})\equiv Li_{2}(\frac{x_{0}}{x_{0}-x_{i}})-Li_{2}(\frac{x_{0}-1}{x_{0}-x_{i}})$, $Li_{2}(z)$ is the di-logarithm function; $x_{1,2}$ are solutions of the equation $x^{2}-\left(\frac{m_{h}^{2}-M_{1}^{2}+M_{2}^{2}}{m_{h}^{2}}\right)x+\frac{M_{2}^{2}-i\delta}{m_{h}^{2}}=0$; $x_{0}=\frac{M_{2}^{2}-M_{0}^{2}}{m_{h}^{2}}$; and $x_{3}=\frac{-M_{0}^{2}+i\delta}{M_{1}^{2}-M_{0}^{2}}$.

Based on denner , the $B$-functions with small absolute values of external momenta can be written in stable forms in numerical computations. Defining $y_{ij}$ ($i,j=1,2$) are solutions of the equation $y^{2}p^{2}-y(p^{2}_{i}+M^{2}_{i}-M^{2}_{0})+M_{i}^{2}-i\delta=0$. New functions $f_{n}(y)$ are defined as follows,

[TABLE]

so that they can be evaluated numerically stable way by choosing

[TABLE]

The $B$-functions now can be expressed in terms of $f_{n}(y)$, namely

[TABLE]

Finally, the $B^{(12)}_{0}$ and $C_{1,2}$ functions are determined as follows,

[TABLE]

In our work above use the following notations, $m_{1}\equiv m_{a}$, $m_{2}\equiv m_{b}$, $p_{1}\equiv p_{a}$ and $p_{2}\equiv p_{b}$.

Appendix B Matching with notations in previous works

This section will show the equivalence given in (23). We recall notations used in EArganda ; iseesaw ; e1612.09290 as follows. The external momenta are $p^{\prime}_{1}$,$(-p_{2}^{\prime})$, and $p_{3}^{\prime}$ for ingoing Higgs boson, outgoing leptons $e_{a}$ and $e_{b}$,respectively. The prime is used to distinguish from the notions that were used in our work, especially those given in Sec. A. Three denominators of the propagators are $D^{\prime}_{0}=k^{2}-m_{1}^{2}$, $D^{\prime}_{1}=(k+p^{\prime}_{2})^{2}-m_{2}^{2}$ and $D^{\prime}_{2}=(k+p^{\prime}_{1}+p^{\prime}_{2})^{2}-m_{3}^{2}$. The one-lopp-three-point functions are defined as,

[TABLE]

The equivalence between above notations with those given in Sec. A are $p^{\prime}_{1}=p_{1}+p_{2}$, $p^{\prime}_{2}=-p_{1}$, $m_{1,2,3}=M_{0,1,2}$. As a result, we get $D^{\prime}_{0,1,2}=D_{0,1,2}$, leading to $C^{\prime}_{0}=C_{0}$ and $C^{\prime}_{\mu}=C_{\mu}$. But the scalar factors $C_{11,12}$ and $C_{1,2}$ are different, namely $C^{\prime}_{\mu}=C_{11}(-p_{1\mu})+C_{12}(p_{1\mu}+p_{2\mu})=(C_{12}-C_{11})p_{1\mu}+C_{12}p_{2\mu}$. Matching this with definition of $C_{\mu}$ defined in Sec. A. We obtain the equivalence for $C_{1,2}$ in (23). Other $B$-functions is proved easily so we omit here.

Appendix C Form factors in unitary gauge for LFVHD

The contribution from diagram in Fig. 1b) to the LFVHD amplitude is

[TABLE]

The final result is

[TABLE]