Radiatively-induced LFV Higgs Decays from Massive ISS Neutrinos

In the inverse seesaw model (ISS), the smallness of the neutrino masses is related to the smallness of a lepton number violating mass term whilst the seesaw scale is naturally close to the TeV scale, which allows for large effects in lepton flavor and universality violating observables. With the ongoing and planned measurements of the Higgs boson properties at the LHC, we found timely to investigate the possibility of having large lepton flavor violating Higgs decay (LFVHD) rates within the context of the ISS, considering the most generic case where three additional pairs of massive right-handed singlet neutrinos are added to the Standard Model particle content. We present a full one-loop computation of the LFVHD rates and analyze in full detail the predictions as functions of the various relevant ISS parameters, which are required to be compatible with the present neutrino data and the present experimental bounds for the three LFV radiative decays, and also consistent with other constraints, like perturbativity of the neutrino Yukawa couplings. At the end, we conclude on the maximum allowed LFVHD rates within the ISS, which may reach maximal values of order $10^{-5}$ for the $H \to e \bar \tau$ and $H \to \mu \bar \tau$ channels, close to the expected future LHC sensitivities.


Introduction
The fact of considering the discovered scalar particle at the CERN-LHC as the Higgs particle of the Standard Model (SM), with very similar properties and a measured mass of m ATLAS h = 125.5 ± 0.6 GeV [1] and m CMS h = 125.7 ± 0.4 GeV [2], has at present reached a broad consensus.
On the other hand, there is also a major consensus about the fact that the SM must be modified to include the neutrino masses and oscillations in agreement with present data [3], which are nowadays quite impressive and call for an explanation from a theoretical framework the ISS context with three extra pairs of RH neutrinos, the Higgs decays into lepton-antilepton pairs H → l klm with k m (LFVHD), which are being currently explored at the LHC [7]. We refer the reader to our main article [8] for more information and details about the full one-loop computation of the LFV partial decay widths, the complete set of references and our full numerical results.

The Inverse Seesaw Model
The ISS supplements the SM with pairs of RH neutrinos, denoted here by ν R and X, with opposite lepton number. We consider a generic model containing three pairs of fermionic singlets, extending the SM Lagrangian with the following neutrino Yukawa interactions and mass terms: where L is the SM lepton doublet, Φ is the SM Higgs doublet, Φ = ıσ 2 Φ * , with σ 2 the corresponding Pauli matrix, Y ν is the 3 × 3 neutrino Yukawa coupling matrix, M R is a lepton number conserving complex 3 × 3 mass matrix, and µ X is a Majorana complex 3 × 3 symmetric mass matrix that violates lepton number conservation by two units. After electroweak symmetry breaking, the 9 × 9 neutrino mass matrix reads, in the electroweak interaction basis (ν C L , ν R , X), with the 3 × 3 Dirac mass matrix given by m D = Y ν Φ , and the Higgs vacuum expectation value taken to be Φ = v = 174 GeV.
In the one generation case, which allows us to illustrate more simply the dependence on the seesaw parameters, there would be just three ISS model parameters, M R , µ X and Y ν , and there would be just three physical eigenstates: one light ν and two heavy N 1 and N 2 . In the limit µ X m D , M R , the mass eigenvalues are given by: with the light neutrino mass m ν being proportional to µ X , thus making it naturally small, and the two heavy masses m N 1 ,N 2 being close to each other. As a consequence in this µ X m D , M R limit, these two nearly degenerate heavy neutrinos combine to form a pseudo-Dirac fermion. A similar pattern of neutrino mass eigenvalues occurs in the three generation case, with one light and two nearly degenerate heavy neutrinos per generation.
We will consider first the simplest ISS scenarios with diagonal µ X and M R matrices (case I). In order to implement easily the compatibility with present neutrino data in these scenarios, we use here the helpful Casas-Ibarra parametrization [9] that can be directly applied to the ISS case, giving: where V is a unitary matrix that diagonalizes M = R is a complex orthogonal matrix that can be written in terms of three arbitrary complex angles θ 1,2,3 . The input ISS parameters that have to be fixed in this case I are the following: m ν 1,2,3 , µ X 1,2,3 , M R 1,2,3 , θ 1,2,3 and the entries of the U PMNS matrix.
On the other hand, given the interesting possibility of decoupling the low energy neutrino physics from the LFV physics in this ISS model by the proper choice of the input parameters, we will take into account specific ISS scenarios with non-diagonal µ X while keeping diagonal M R (case II) which can provide the largest LFV Higgs decay rates. Once some specific inputs are provided for Y ν and M R , the proper µ X matrix that ensures the agreement between low energy neutrino predictions and data can be easily obtained by solving ISS equations, which leads to: with m D = vY ν and m ν = diag(m ν 1 , m ν 2 , m ν 3 ). It should be noted that for a generic Y ν texture, this µ X will be in general non-diagonal. Therefore, the most relevant input ISS parameters in the case II are Y ν and M R (in addition to the input parameters that are relevant for low energy neutrino physics, m ν and U PMNS ).

Remarks on the LFVHD Computation
In our one-loop computation of the LFV rates we work in the mass basis for all the particles involved, with diagonal charged leptons, and take into account the contributions from all the nine physical neutrinos. As for the gauge choice, we choose the Feynmant'Hooft gauge. The full set of contributing one-loop diagrams can be found in [8], and we have adapted the complete one-loop formulas for the Γ(H → l klm ) partial decay width, taken from [10], to the ISS case. We have focused on the decays H → µτ, eτ, eμ and have not considered their related CP conjugate decays H → τμ, τē, µē which, in the presence of complex phases, could lead to different rates. All the formulas for the LFVHD have been implemented into our private Mathematica code. For the numerical predictions of the BR(H → l klm ) rates, we use m H = 126 GeV and its corresponding SM total width is computed with FeynHiggs [11,12,13] including two-loop corrections.
In order to ensure the validity of the Casas-Ibarra parametrization, we have imposed that the error on the light neutrino masses estimated with it, meaning the differences between the input m ν 1,2,3 and the output m n 1,2,3 masses, was below 10% and that the 9 × 9 rotation matrix exhibited the required unitarity property. Furthermore, since a given set of input parameters can generate arbitrarily large Yukawa couplings, we enforce their perturbativity by setting an upper limit on the entries of the neutrino Yukawa coupling matrix, given by for i, j = 1, 2, 3. At the same time that we analyze the LFVHD, we also compute the one-loop l m → l k γ decay rates within this same ISS framework (using the analytical formulas provided in [14] and [15]) and for the same input parameters, and check that these radiative decay rates are compatible with their present experimental 90% CL upper bounds: Regarding the Higgs total width, we focus on the scenario where the new fermionic singlets have a mass above 200 GeV, thus they do not contribute to new invisible decays. To avoid potential constraints from lepton electric dipole moments, we assume that all mass matrices are real, as well as the PMNS matrix. Additional constraints might also arise from lepton universality tests [18,19]. Nevertheless, in the scenario that we consider where the sterile neutrinos are heavier than the Higgs boson, points that would be excluded by lepton universality tests are already excluded by µ → eγ. At the end, we found that the most stringent constraints for our study are by far µ → eγ and the Yukawa coupling perturbativity limit of Eq. (7).

LFV Rates in the Inverse Seesaw: Case I
We study first the LFV rates as functions of the most relevant ISS parameters within the case I, trying to localize the areas of the parameter space where the LFVHD can be both large and respect the constraints on the LFV radiative decays. The results will be presented in two generically different scenarios for the heavy neutrinos: 1) the case of (nearly) degenerate heavy neutrinos, and 2) the case of hierarchical heavy neutrinos.
The case of (nearly) degenerate heavy neutrinos is implemented here by choosing degenerate entries in M R i = M R and µ X i = µ X (i = 1, 2, 3). First we show in Fig. 1 the results for all LFV rates as functions of the common RH neutrino mass parameter M R for all LFVHD channels (upper panel) and for the LFV radiative decay channels (lower panel). As expected, we find that the largest LFVHD rates are for BR(H → µτ) and the largest radiative decay rates are for BR(τ → µγ).
We also see that, for this particular choice of input parameters, all the predictions for the LFVHD are allowed by the present experimental upper bounds on the three radiative decays for all explored values of M R in the interval (200, 10 7 ) GeV. Besides, it shows clearly that the most constraining radiative decay at present is by far the µ → eγ decay. This is so in all the cases explored in this work, so whenever we wish to conclude on the allowed LFVHD rates we will focus mainly on this radiative channel.
Regarding the M R dependence shown in Fig. 1, we clearly see that the LFVHD rates grow faster with M R than the radiative decays which indeed tend to a constant value for M R above ∼ 10 3 GeV. In fact, the LFVHD rates can reach quite sizable values in the large M R region of these plots, and yet stay allowed by the constraints on the radiative decays. M R of the LFVHD and the radiative rates shown by Fig. 1 is an interesting feature that we wish to explore further. As it is clearly illustrated in Fig. 2, the radiative decay rates can be well approximated for large M R by a simple function of |(Y ν Y † ν ) km | 2 given by: which provides very close predictions to the exact rates for M R > 10 3 GeV. Then we can understand the final constant behavior of all the radiative decay rates with M R , since the |(Y ν Y † ν ) km | 2 elements grow with M R approximately as M 4 R in the parametrization used here of Eq. (5). This simple behavior with M R is certainly not the case of the LFVHD rates, and we conclude that these do not follow this same behavior with |(Y ν Y † ν ) km | 2 . In order to understand the different functional behavior of LFVHD rates with M R , we have tried to find an approximate simple formula that could explain the main features of these rates. As we have already said, in contrast to what we have seen for the LFV radiative decays in Eq. (11), a simple functional dependence being proportional to |(Y ν Y † ν ) 23 | 2 is not enough to describe our results for the BR(H → µτ) rates. Considering that, in the region where the Yukawa couplings are large, we have looked for a simple expression that could fit properly the contributions from the dominant diagrams. From this fit we have found the following approximate formula: which turns out to work reasonably well. In Fig. 3 we show the predicted rates of BR(H → µτ) with 1) the full one-loop formulas; 2) taking just the contributions from dominant diagrams, and 3) using Eq. (12) (dotted lines). We see clearly that this Eq. (12) reproduces extremely well the contributions from dominant diagrams and approximates reasonably well the full rates. The approximation is pretty good indeed for the M R region above the dips. The change of functional behavior with M R in the two different M R regions, from nearly flat with M R in the approximate result to fast growing as ∼ M 4 R , also gives a reasonable approach to the full result, as well as the appearance of dips. The location of the dips is however not so accurately described by the approximate formula, since in the region where the cancellation among the dominant diagrams takes place the other diagrams (not considered in the fit) also contribute. Overall we find the approximate formula given by Eq. (12) very useful for generic estimates in the ISS, which could be also applied to other parametrizations of the neutrino Yukawa couplings.
Next we study the dependence of the LFV rates on µ X . The behavior of BR(H → µτ) and BR(µ → eγ) as functions of µ X , for several values of M R , m ν 1 = 0.1 eV and R = I, are displayed in Fig. 4. Both LFV rates decrease as µ X grows; however, the functional dependence is not the same. The LFV radiative decay rates decrease as µ −2 X , in agreement with the approximate expression (11), while the LFVHD rates go as µ −4 X when the Yukawa couplings are large. For a fixed value of µ X , the larger M R is, the larger BR(H → µτ) can be reached, while the same prediction for BR(µ → eγ) is obtained for any value of M R . We have already learnt this independence of the LFV radiative decays on M R from the previous figure, which can be easily confirmed on the lower panel of Fig. 4, where all the lines for different values of M R are superimposed. We observe again the existence of dips in the upper panel of Fig. 4. We see in this figure that the smallest value of µ X allowed by the BR(µ → eγ) upper bound is µ X ∼ 5 × 10 −8 GeV, which is directly translated to a maximum allowed value of BR(H → µτ) ∼ 10 −11 , for M R = 10 4 GeV.
Once we have studied the behavior of all the LFV observables considered here with the most relevant parameters, we next present the results for the maximum allowed LFVHD rates in the case of heavy degenerate neutrinos. The plot in Fig. 5 shows the contour lines of BR(H → µτ) in the (M R , µ X ) plane for R = I and m ν 1 = 0.1 eV. These contour lines summarize the previously learnt behavior with M R and µ X , which lead to the largest values for the LFVHD rates in the bottom right-hand corner of the plot, i.e. at large M R and small µ X . We also notice the appearance of dips in the (M R , µ X ) plane which correspond to the dips commented before in the previous figures. The most important conclusion from this contour plot is that the maximum allowed LFVHD rate for this simple hypothesis of diagonal and degenerate µ X and M R matrices is approximately BR(H → µτ) ∼ 10 −10 and it is found for M R ∼ 2 × 10 4 GeV and µ X ∼ 5 × 10 −8 GeV. We have found similar conclusions for BR(H → eτ). The case of hierarchical heavy neutrinos refers here to hierarchical masses among generations and it is implemented by choosing hierarchical entries in the M R = diag(M R 1 , M R 2 , M R 3 ) matrix. As for the µ X = diag(µ X 1 , µ X 2 , µ X 3 ) matrix which introduces the tiny splitting within the heavy masses in the same generation we choose it here to be degenerate, µ X 1,2,3 = µ X . We focus on the normal hierarchy M R 1 < M R 2 < M R 3 , since we have found similar conclusions for other hierarchies. Figure 6 shows the contour lines of BR(H → µτ) in the (M R 3 , µ X ) plane for R = I and m ν 1 = 0.1 eV. It is clear from this contour plot that the behavior of the LFV rates in the hierarchical case with respect to the heaviest neutrino mass M R 3 is very similar to the one found previously for the degenerate case with respect the common M R . Again, there are dips in the BR(H → µτ) rates due to the destructive interferences among the contributing diagrams. We have found, for this hypothesis of degenerate µ X and hierarchical M R , an enhancement of the LFVHD rates by approximately one order of magnitude as compared to the degenerate case in most of the explored parameter space regions, where the maximum allowed BR(H → µτ) rates reach values up to about 10 −9 for M R 1 = 900 GeV, M R 2 = 1000 GeV, M R 3 = 3 × 10 4 GeV, µ X = 10 −7 GeV, and R = I.

LFV Rates in the Inverse Seesaw: Case II
In the line of thinking of increasing H → eτ and H → µτ rates, while reducing the predictions of µ → eγ rates, which is the most restrictive constraint together with the Yukawa perturbativity, we propose in this section more general scenarios with non-diagonal µ X and diagonal and degenerate M R . In order to localize the class of scenarios leading to large and allowed LFVHD rates, we first make a rough estimate of the expected maximal rates for the H → µτ channel by using our approximate formula of Eq. (12) which is given just in terms of the neutrino Yukawa coupling matrix Y ν and M R . The most constraining observables in our study are the LFV radiative decays, which limit the maximal value of the non-diagonal (Y ν Y † ν ) i j entries. By using our approximate formula of Eq. (11) and the present bounds in Eq. (10), we obtain: Then, in order to simplify our search, and given the above relative strong suppression of the 12 element, it seems reasonable to neglect it against the other offdiagonal elements. In that case, Eq. (12) for the H → µτ decay mode reads as BR approx H→µτ This equation clearly shows that the maximal BR(H → µτ) rates are obtained for the maximum allowed values of (Y ν Y † ν ) 23 , (Y ν Y † ν ) 22 and (Y ν Y † ν ) 33 . By setting the maximum allowed value for v 2 (Y ν Y † ν ) max 23 /M 2 R to that given in Eq. (15) and fixing the values of (Y ν Y † ν ) 22 and (Y ν Y † ν ) 33 to their maximum values allowed by Eq. (7), we obtain our approximate prediction for the maximal rates: BR max The same procedure can be applied to the H → eτ channel, leading to similar results on its maximal branching ratios. These maximal H → µτ and H → eτ rates are much more promising than the ones obtained in the previous section under the simple hypothesis of diagonal µ X and M R , since they are closer to the expected LHC sensitivity for these channels (see for instance [20,21]) which will be surely improved in the forthcoming runs of the LHC. Assuming this kind of Yukawa matrix which gives rise to the maximal BR(H → µτ) and BR(H → eτ), the issue now would be to find the possible µ X matrices that keep the agreement with neutrino data. This fact will be always possible by means of Eq. (6) and because of the possibility of decoupling the low energy neutrino physics from the LFV physics in the ISS model. This means that, for a given input Y ν , Eq. (6) tells us which µ X keeps the agreement with neutrino data. We address the reader to our main article [8] for more details on this case II and for some illustrative examples of Yukawa textures producing these maximal LFVHD rates.

Conclusions
We have studied the LFVHD within the context of the ISS where three additional pairs of massive RH singlet neutrinos are added to the SM particle content. The most relevant ISS parameters have been found to be M R and Y ν . We have required that the input ISS parameters be compatible with the present neutrino data and other constraints, with the perturbativity of the neutrino Yukawa couplings and the present upper bounds on the LFV radiative decays being the most restrictive ones. The largest maximum LFVHD rates are BR(H → eτ) and BR(H → µτ) and under the hypothesis of diagonal µ X and M R they can reach at the most values of 10 −10 for the degenerate heavy neutrino case and 10 −9 for the hierarchical case. Going beyond this simple hypothesis, we find more general ISS scenarios with non-diagonal µ X in which BR(µ → eγ) is extremely suppressed and BR(H → eτ) or BR(H → µτ) are larger while respecting the BR(τ → eγ) and BR(τ → µγ) upper bounds. These rates reach maximal values of 10 −5 , being very promising for LFVHD searches at the LHC.