Heavy neutral leptons — Advancing into the PeV domain

: Heavy neutral leptons (HNLs) are hypothetical particles able to explain neutrino oscillations and provide a mechanism for generating the baryon asymmetry of the Universe. Quantum corrections due to such particles give rise to ﬂavor violating processes in the charged lepton sector. Based on the fact that these corrections grow with HNL masses, we improve existing constraints by orders of magnitude in mass and mixing angle. This allows us to probe part of the parameter space of leptogenesis with multi-TeV HNLs. We also show that one will be able to infer HNL parameters in a signiﬁcant portion of the parameter space for TeV-PeV masses if charged lepton ﬂavor violating signals are detected.


Introduction
Neutrinos in the Standard Model come in three different flavors.However, unlike other fermions, neutrinos change their flavor when propagating -a phenomenon known as neutrino oscillations (see e.g.[1,Ch. 14]).The Standard Model forbids such a flavor-changing transition.Therefore neutrino oscillations imply the existence of extra quantum states, see e.g.[2] for an overview of neutrino mass models.One example of such states is given by right-chiral gauge singlet fermions, interacting with lepton doublets and the Higgs field via the Yukawa interactions.Being gauge singlets, these states can have Majorana masses, whose scale can be arbitrary [3][4][5][6][7][8][9][10]. Compared to the Standard Model spectrum, this Type I seesaw model contains light neutrinos, ν i , i = 1, 2, 3 whose masses and mixings fit the experimental neutrino measurements [11] and, additionally, heavy neutral leptons (HNLs) N I , I = 1, 2 . . .that interact similarly to neutrinos, but are much heavier and have a suppressed interaction strength due to the mixing angle U2 α . 1 HNL masses can lie anywhere from O(1 eV) till O(10 15 GeV), with experiments limiting their interaction from above (direct experimental searches, see [12,13]) or from below (as in the case of primordial nucleosynthesis, see e.g.[14][15][16]).
In this work, we highlight the non-trivial mass dependence of the cLFV processes and demonstrate that HNLs within the TeV-PeV mass range could significantly influence the magnitude and rates of cLFV events.More specifically, we show that for HNLs with masses M N ≫ 1 TeV, the existing cLFV bounds are much stronger than previously estimated.In fact, they are one of the the strongest across the whole range of HNL masses. 2  This observation suggests the feasibility of using precision frontier measurements as a strategy for investigating the elusive high-mass domain of HNLs.We show that it is possible to explore a portion of the leptogenesis parameter space previously considered inaccessible.Furthermore, we demonstrate that potential observations of cLFV signals will allow for a more accurate determination of HNL parameters.

Basic definitions for heavy neutral leptons
In this section, we collect the basic definitions related to the Type-I Seesaw Lagrangian.While it contains only known results, they are needed for the consistency of our exposition.

Type I seesaw Lagrangian
The Lagrangian of the type-I seesaw model reads: where L SM is the SM Lagrangian, ν RI , I = 1, . . ., N are the right-chiral components of neutrino fields -SM gauge singlets -, L α is SM lepton doublets, Φ = iσ 2 Φ where Φ is the Higgs doublet, F αI is a Yukawa matrix, and M I is the Majorana masses of righthanded neutrinos.Finally, sum over I and α is assumed in (2.1).After spontaneous symmetry breaking, the model will acquire Dirac masses, (m D ) αI mixing right and lefthanded neutrino components.In the case |m D | ≪ M I , the diagonalization of the mass matrix results in three fields ν i with light masses, m i and N heavy fields N I with heavy masses M I of the shape: where V pmns is the PMNS matrix [74], and α = e, µ, τ .Eq. (2.2) is the famous seesaw mechanism [3,4,[7][8][9][10].
After diagonalization, neutrinos in the flavor basis are in a superposition of light and heavy mass states: where Θ is the left-right neutrino mixing angle, and the 3 × (3 + N ) matrix U is defined as Eq. (2.4) makes it clear that HNLs will interact in the exact same way as active neutrinos do, with the notable exception that it will be suppressed by the mixing angle Θ αI forming a 3 × N matrix with components Furthermore, we can define flavor mixing angles to which in the case of nearly degenerate HNLs the signal would be proportional.
In addition, after diagonalization interactions between neutral leptons and SM gauge bosons get modified from the usual SM ones.These interactions read as: where We properly used the Feynman rules regarding Majorana particles [75], especially with regards to vertices involving lepton number violating (LNV) processes or in interactions with neutral bosons that involve two Majorana particles.
There is a number of relations that the U, C and mass matrix of neutral fermions follow that can be found in Appendix A.
The model has introduced a total of new 7N − 3 different parameters to the SM.Neutrino oscillation data already constraints some of them: two different mass splittings and three mixing angles from the PMNS matrix [11].In the active neutrino sector, we additionally have the lightest neutrino mass, two Majorana CP phases, δ CP , one Dirac CP phase, and the discrete choice of mass hierarchy.
Throughout this work, we are dealing with the simplest realistic case of N = 2.In this case, the lightest neutrino mass is set automatically to zero, and only a linear combination of two Majorana phases contributes, which we will call η.This leaves us with three unknown parameters in the active neutrino sector, with the already mentioned unbeknownst mass hierarchy.In this model, the R matrix depends on the choice of hierarchy [77,78]: It is useful to define the sum of the absolute value squared of all the elements of |Θ αI | 2 , U 2 tot , in terms of Casas-Ibarra parameters, and HNL and active neutrino masses.The expression becomes simple when considering two HNLs with equal masses: the reader may find a more general expression in [79] for non-degenerate HNL masses.
3 Decoupling behaviour cLFV observables receive contributions from HNLs, appearing in the loop once or twice: compare the diagram (a) with the diagrams (b) and (c) in Figure 1 (see Appendix B for details).Two-HNL contributions bear similarity with the "penguin diagrams" in flavor physics [see e.g., 80] and appear whenever Z or Higgs bosons contribute to the diagrams.Schematically, the conversion rate µ → e on a nucleus (the process that will give the strongest bounds) has the form where The functions c i (M N ) scale for HNL mass M N ≫ m W as: As a result, although subdominant due to U 2 tot ≪ 1, the strong mass dependence of c 2 (M N ) makes the second term larger than the first one.As a result, HNL contributions to cLFV processes become much stronger at large masses than previously estimated. 3he above result seemingly contradicts the "Decoupling Theorem" [81], as loop effects of heavy particles grow with mass.However, theories that undergo spontaneous symmetry breaking violate this theorem if dimensionless couplings -in our case, the Yukawa coupling that determines the Dirac mass -are sent to infinity (see e.g.[82][83][84][85][86], or the book [87,Chapter 8]).At a technical level, this happens because the coupling between HNLs and Goldstone bosons is proportional to the Yukawa couplings, and in turn, to the Dirac Mass term m D ∼ ΘM N (see Eq. 2.6).If one considers Θ and M N as independent parameters, then Dirac mass term grows as M N → ∞, and there is no decoupling.On the other hand, keeping Yukawa constants fixed implies |Θ| ∝ M −1 N .As a result, the effect (3.1) disappears.By increasing M N at fixed |Θ| the HNL decay width Γ N ∝ |Θ| 2 M 3 N exceeds 1 2 M N (equivalently, Yukawa coupling constant of neutrino exceeds √ 4π).The corresponding non-perturbative region is shaded in gray below.See [88][89][90][91], and Appendix D for a brief discussion on the stability on second-order terms.Indeed, the non-decoupling behaviour of HNLs implies that when we match the amplitudes to a series of EFT operators, we will find that such operators will depend on HNL masses or on their Yukawa couplings, and this is exactly what happens when one does a one-loop matching to different EFT operators [92].

Results
Using the proper expression for the cLFV branching ratios/conversion rates, we reinterpret the existing limits on HNL parameters as shown in Figure 2 (solid lines).To efficiently present and compare our results, we fix ratios of U 2 α to the benchmark values, dictated by neutrino oscillations [12]: 0.06 : 0.48 : 0.46 (NO) 0.33 : 0.33 : 0.33 (IO) (4.1) (see [96]) and express all the limits in terms of the total mixing angle, U 2 tot .Compared to the electroweak precision data constraints [69], our bounds improve by as much as three orders of magnitude.Compared to the previous cLFV constrains (not taking into account the c 2 (M N ) terms), we improve the bounds by roughly one order of magnitude for masses O(100 TeV).To the best of our knowledge, the constraints from µ → e conversion in gold [93] have not been previously used to derive bounds on HNL parameters in the multi-TeV region, and the only work that took into account c 2 (M N ) terms when deriving the cLFV limits was [68] that however did not explore them beyond ∼ 10 TeV.
The constraints from µ → 3e and muon conversion in the case of HNLs with a hierarchical mass spectrum are presented Appendix B.

Indirect detection of TeV-PeV HNLs
In the wide range of HNL parameters, two or more future cLFV experiments can detect a signal, as Figure 4 demonstrates.Non-trivial mass dependences of µ → e conversion rates and µ → 3e processes open a possibility to recover HNL parameters in this case.Similar idea has been explored in the past, [see e.g.94, 105, 106], albeit without "penguin" contributions that change the mass dependence.
Consider the simplest HNL model explaining neutrino oscillation data -two HNLs with almost degenerate masses. 4The cLFV branching ratios (3.1) depend in this case on three parameters: common mass, M N , total mixing angle, U 2 tot , and Λ eµ .While flavor mixing angles, U 2 α can still vary in the wide range, Λ eµ changes only weakly if the mass ordering and the CP phase is known [see e.g.111,112].As a result, given two cLFV measurements (µ → e conversions and µ → 3e measurement, Table 1), the HNL parameters can be reconstructed.M N (GeV)  , see text for details.The black line in the right plot is the upper boundary of the parameter space where successful leptogenesis with 3 HNLs is possible [28] (dashed-dotted part is an extrapolation to guide the eye).The sensitivity reach of several future experiments (SHiP [103], FCC-ee [104], FCC-hh [48]) is shown as thin dashed lines.To demonstrate the feasibility of parameter reconstruction, we have randomly assumed HNL parameters in the light green region of Figure 4 and computed the rates of the cLFV processes.Assuming that these were the measured rates (with 50% measurement error), we defined the log-likelihood as a mean square deviation between the predicted and "measured" rates.We performed an MCMC scan of the parameter space using the emcee pacakge [113].Using Kernel Density Estimation method implemented in GetDist [114], we successfully recovered both mass and the total mixing angle, see Figure 5. .Reconstruction of HNL parameters in case of two cLFV detections.HNL parameters, explaining neutrino data [11], assuming normal ordering, and that the provided cLFV signals are within the future discovery zone (c.f. Figure 4).We built a likelihood function, depending on four parameters (shown in the plots as straight lines) and explored it with MCMC sampling.ω is the mixing angle from the Casas-Ibarra matrix and η is the Majorana Phase (see Sec. 2.2).Dotted grey lines show the "true" parameters.

Discussion
Heavy neutral leptons couple to neutrinos of different flavors and thus give rise to charged lepton flavor violating processes.Many laboratories worldwide (Table 1) pursue searches for these cLFV processes.Their negative results provide constraints on HNL parameters, see e.g.[69].For µ → e conversion rate, the HNL contributions to cLFV processes contain terms proportional to the U 4 tot and those, proportional to the U 8 tot (with the common factor Λ eµ , Eq. (3.2)).The latter terms have been neglected in the past.However, for M N ≫ m W , the term U 8 tot is multiplied by the growing function in mass, (M N /m W ) 4 .Therefore, for M N ≫ m W this second term becomes dominant, making constraints stronger by orders of magnitude, as we showed in this work (see Figure 2 for current and Figure 3 for future sensitivity).The resulting constraints for realistic HNL models (i.e.those where HNLs mediate neutrino oscillations) allow us to reach limits on the mixing angle compared to Table 1.Present cLFV limits and future projections for various Λ αβ processes.
the best direct searches for masses up to ∼ 30 PeV.In particular, these constraints exceed by orders of magnitude in mass reach and coupling sensitivity those that can be achieved with direct searches for HNLs at FCC-hh [13,[47][48][49].We derived our main results for the case of the 2 HNL model, but they are readily generalizable to the models with 3 HNLs (see Appendix B).In this case, cLFV measurements can probe part of the parameter space of the leptogenesis [28] as Figure 3 demonstrates.See also [134].
The present analysis has assumed benchmark flavor ratios (4.1).Scanning over the whole allowed range of Λ αβ does not qualitatively change our results (see also [135]).In particular, in models with 2 HNLs, even in the case of normal ordering, where U 2 e ≪ U 2 µ,τ , the processes involving LFV in the eµ sector dominate the current and future limits.In the models with 3 HNLs the suppression of the electron mixing may become even more extreme [13], potentially making µτ channels competitive.This demonstrates importance of the future searches, such as TauFV [136][137][138] at the CERN Beam Dump Facility [139] measurement τ → 3µ processes at the level ∼ 10 −10 .

A Unitarity and Seesaw relations
The matrix elements U αi of the matrix U, defined by Eq. (2.5), obey the following unitarity and seesaw relations [56]: Relation (A.1) is the unitarity of the U matrix, whereas (A.2) defines the coupling between the Z boson and n i n j .The equalities (A.3) are a consequence of two previous ones, and the last set of equivalences are nothing but the Seesaw relation, expressed in Eq. (2.2), but written in terms of U and C.
To include the "HNL penguin" diagrams, we modified the FeynRules model file HeavyN [47,48,150] and [151] and added previously missing vertices involving the coupling between two HNLs and neutral bosons -those in Eq. (2.9), (2.11) and (2.10).The updated model is available at [152].3).Both parts are asymptotic in mass, the two of them get two powers in mass from the couplings between HNLs and Goldstone bosons.The former also gets two additional powers in mass coming from the propagators of the HNLs, which cancels the mass suppression coming from the loop functions.The latter gets the momentum from the HNL propagators, but the leading terms of the loop functions are not suppressed in mass. 5Overall, both amplitudes will have a mass squared dependence.
The origin of the asymptotic behavior of diagrams (e) and (i) in Fig. (C.9) can be explain similarly.Both diagrams get four powers in mass from the couplings between the HNLs and Goldstone bosons.Diagram (e) picks up the momentum terms from the propagators of the HNLs, this gives it a loop function which removes two powers of mass from the amplitude.Whereas, diagram (i) picks up the mass terms from the propagators, but has a stronger suppression from the loop functions.All in all, both terms will have a mass squared dependence.
The diagrams that have a red letter attached to them are divergent.All of these divergences vanish when we sum up the diagrams.Namely, diagrams (c, d, e) in Figures C. 7 It should also be noted that in the limit when the four-momentum squared of the outer boson is zero, Q 2 → 0, the sum of the three diagrams vanishes completely.The reason for this becomes clear when the calculation is done in the R ξ gauge, where the three diagrams are completely gauge-dependent.

B.2 Branching ratio for cLFV processes
The cLFV processes that receive contributions from "HNL penguins" or "HNL box diagrams" include ℓ β → ℓ α ℓ γ lη .The formulas of the branching ratios are different depending on whether α = γ = η, α = η ̸ = γ, or α = γ ̸ = η.In the first case, there is an additional factor coming from identical fermions in the final state that is not found in the second case; and in the third case, only box diagrams contribute, Z and γ penguins do not contribute at 1-loop.
The formulas of the branching ratios for these decays are: Another process of interest, µ → e conversion on the nucleus A, also receives "penguin" contributions.The conversion rate can be computed using the results obtained in [153].The conversion rate is: (B.4) we also mention for completeness that the processes ℓ β → ℓ α γ do not contain "HNL penguins" and therefore do not exhibit the behavior in question: The loop functions that enter expressions (B.1-B.5) are given by where , where g is the weak coupling constant, e = √ 4πα is the value of the elementary charge, and A , and are all constants which depend on the nucleus in question, their precise definitions and values for difference nuclei can be found in [153].
Equations (B.6-B.11)can be rewritten in terms of the light-heavy mixing angle, Θ and HNL masses alone, using the fact that U U † = 1 and that neutrino masses are completely negligible in the limits we're considering: The functions G γ , F γ , F Z , F Xbox , F box and G box are listed in the appendix B.4 below.Finally, another set of decays that also receives contributions from HNL penguins are Z → ℓ α lβ and H → ℓ α lβ [56-58, 62, 154-156], but neither of their constraints are competitive with the other processes mentioned.Moreover, Higgs decays to two leptons are helicity suppressed: their branching ratio is proportional to the masses of the outgoing particles, adding another suppression to it.

B.3 Relevant combinations of the mixing angles
We see that the relevant observables depend on specific combinations of products of the elements of Θ.Specifically, the branching ratios are proportional to the parameter Λ αβ , defined by Eq. 2 in the main text).HNL penguins depend on the combinations In the case of two HNLs with equal masses and Im(ω) ≫ 1, this expression simplifies to It should be noted that there is another kind of HNL penguin diagrams, where each HNL, running in the loop, violates the total lepton number by +1 or −1.The resulting process is still lepton-number conserving but is proportional to the sum (B.18) albeit with C * IJ instead of C IJ .In this case, the sum automatically vanishes.
Another noteworthy combination involves a sum that uses the same flavor thrice, relevant for ℓ α → 3ℓ β decays.It can be easily evaluated using the same limits as we used to evaluate the last one: where λ β = U 2 β /U 2 tot , with U 2 β defined in Eq. (2.7).Moreover, considering the same case, the sum I,J Θ αI Θ * βJ Θ βI Θ * βJ would immediately vanish.These results are valid for any hierarchy of HNL masses.But it can be seen from Eq. (B.14) and Eq.(B.15) that the sums are also proportional to HNL masses.In the case of degenerate HNL masses, we can take all the mass-dependent parts out of the summation, and we would only have to deal with the sum of mixing angles.If HNL masses were not degenerate, then the LNV parts would come to contribute. Figure (B.6)shows how much the constraints would change if we had hierarchical HNLs.
At this point, we have to stress that including hierarchical HNLs no longer protects the smallness of active neutrino masses, as there would no longer be a lepton number conserving symmetry guaranteeing their smallness.Indeed, if HNLs were to be highly hierarchical, loop contributions to neutrino masses would become larger compared to experimental results, rendering the theory unnatural [32,34,35,157,158].

B.4 Loop functions
For completeness we also show the relevant loop integrals that appear in the expression (B.12-B.17): We have only considered the leading terms in x for each power of U tot in equations (B.28-B.29),there are additional powers of x that are not negligible and should be taken into consideration when doing any calculation regarding them.Equations (B.2-B.3) also show a similar asymptotic behaviour to that of (B.1), as they also have contributions coming from Z-penguins and HNL box diagrams.The decays these formulas describe (τ → µµe, τ → eeµ) are not as constrained as muon decays and were not examined, but they do offer an interesting opportunity to probe all mixing angles since these decays require HNLs to couple to all lepton flavors.

C Feynman diagrams contributing to the relevant processes
We show below all the relevant Feynman diagrams of the processes we highlighted.All are in the Feynman-'t Hooft gauge and therefore contain intermediate Goldstone bosons.The diagrams containing divergences have their respective letter in red, and the ones that exhibit the non-decoupling behavior highlighted in the main text are in green.where the 16π 2 is the usual term one gets when doing the four dimensional loop integrals, which in case of a two-loop diagram it should be squared, of course.This is in direct comparison with the parametric dependence of the leading order terms of one-loop diagrams: The two-loop amplitude will be smaller than the one-loop one as long as Y < 4π, which is in complete accordance with the classical perturbative limit for Yukawa couplings.One should expect the higher order ones to follow a similar dependence on Yukawas and a suppression of 16π 2 .

Figure 1 .
Figure 1.One-loop diagrams contributing to the µ → e conversion process, including 2 HNL "penguin" diagrams (a).Diagrams containing HNLs that dress external lines are not shown (see Appendix C for the full list of such diagrams).

Figure 2 .
Figure2.Improved experimental limits on HNLs coupling from charged lepton flavor violation processes (solid lines) for the normal ordering (left) and inverted ordering (right).The strongest bound comes from the µ → e conversion in gold[93].Blue shaded region: previous constraints, not including cLFV.Gray shaded region: non-perturbative regime (see text for details).Dotted lines: previous cLFV constrains from the works[65,68,69,94,95].Results are expressed in terms of U 2 tot

Figure 3 .
Figure 3. Constrains on HNL parameters from future cLFV experiments for the normal Ordering (left) and inverted Ordering (right).Two HNL contributions are taken into account where applicable.Notice the different range of y-axes as compared to Figure 2. Maximal probed mass may reach 3 × 10 7 GeV should the recently proposed µ → e conversion experiment [102] be performed.Results are expressed in terms of U 2 tot for flavor ratios (4.1), see text for details.The black line in the right plot is the upper boundary of the parameter space where successful leptogenesis with 3 HNLs is possible[28] (dashed-dotted part is an extrapolation to guide the eye).The sensitivity reach of several future experiments (SHiP[103], FCC-ee[104], FCC-hh[48]) is shown as thin dashed lines.

Figure 4 .
Figure 4. Parameter space probed by future cLFV experiments for the normal ordering (left) and inverted ordering (right).Color bands illustrate where one or several future cLFV experiments have sensitivity and thus can observe a signal.Observation of a signal in two or more experiments (orange to green bands) will allow for reconstructing parameters of realistic HNL models.The white region is explored by the current cLFV experiments.

Figure 5
Figure 5. Reconstruction of HNL parameters in case of two cLFV detections.HNL parameters, explaining neutrino data[11], assuming normal ordering, and that the provided cLFV signals are within the future discovery zone (c.f.Figure4).We built a likelihood function, depending on four parameters (shown in the plots as straight lines) and explored it with MCMC sampling.ω is the mixing angle from the Casas-Ibarra matrix and η is the Majorana Phase (see Sec. 2.2).Dotted grey lines show the "true" parameters.
The diagrams that display the apropos asymptotic behavior in mass are shown below in Appendix C in green.The anomalous Zℓ α ℓ β coupling shown in Fig. (C.8) also gives a non-decoupling effect on both ℓ β → 3ℓ α and muon conversion in a nucleus processes.The amplitude of diagram (j) in Fig. (C.8) contains two different parts, one dependent on C ij and another on C * ij , the latter is suppressed in the case of degenerate HNLs (see B. -C.8 have their divergences vanish due to the unitarity relation expressed in Eq. (A.1), this is also the case for Diagram (i) in Fig. C.8. Whereas for Diagram (j) in Fig. (C.8 vanishes due to the Seesaw relation, expressed in Eq. (A.4).Diagrams (f, g, h) in both, Figure C.7 and Figure C.8 vanish when the three diagrams are summed.

Figure B. 6 .
Figure B.6.Constraints for a degenerate and hierarchical HNL mass spectrum for two different processes (showed at the top of the plot): muon conversion in Gold (left) and µ → 3e (right).

B. 5
Asymptotic behaviour of branching ratiosUsing the formulas defined in Appendix B.4, we can extract the behaviour of branching ratios (B.1-B.4) for large masses x

Figure C. 7 .
Figure C.7. Diagrams contributing to the γℓ α ℓ β vertex.To this vertex there are no two-HNL contributions.The diagrams with red labels have an UV divergence that ultimately vanishes.

Figure C. 8 .Figure C. 9 .Figure C. 10 .Figure D. 11 .16π 2 2 ,
Figure C.8. Diagrams contributing to the Zℓ α ℓ β vertex.The diagrams in green are new, accounting for two-HNL interactions.The diagrams with red labels have an UV divergence that ultimately vanishes.