Effective portals to heavy neutral leptons

The existence of right-handed neutrinos, or heavy neutral leptons (HNLs), is strongly motivated by the observation of neutrino masses and mixing. The mass of these new particles could lie below the electroweak scale, making them accessible to low-energy laboratory experiments. Additional new physics at high energies can mediate new interactions between the Standard Model particles and HNLs, and is most conveniently parametrized by the neutrino Standard Model Effective Field Theory, or $\nu$SMEFT for short. In this work, we consider the dimension six $\nu$SMEFT operators involving one HNL field in the mass range of $\mathcal{O}(1)$ MeV $<M_N<\mathcal{O}(100)$ GeV. By recasting existing experimental limits on the production and decay of new light particles, we constrain the Wilson coefficients and new physics scale of each operator as a function of the HNL mass.


Introduction
The evidence for neutrino masses and mixing, as observed in neutrino oscillations, calls for an extension of the Standard Model (SM) of particle physics. Among the many different options, the completion of the SM particle content with the right-handed counterpart of the left-handed neutrinos is the simplest possibility. Apart from the usual Yukawa couplings with the Higgs, the unique gauge-singlet nature of right-handed neutrinos allows for them also a lepton-number-violating Majorana mass term. This new scale would be the first dimensionful term in the Lagrangian not directly related to the electroweak scale and the Higgs mechanism, and, as such, its value remains utterly unknown.
A common assumption is that the Majorana mass of the right-handed neutrinos is much larger than the electroweak scale. This heavy scale would suppress the induced mass for the mostly-active neutrino mass eigenstates via the d = 5 Weinberg operator [1], providing a rationale for their smallness. This is the celebrated high-scale Seesaw mechanism [2][3][4][5] which, as a bonus, may also explain the origin of the matter-antimatter asymmetry of the universe via the leptogenesis mechanism [6]. The drawback of this scenario is that the same very large mass scale that suppresses the light neutrino masses also renders the model virtually untestable since all phenomenology would be similarly suppressed. Furthermore, the addition of a new mass scale much above the electroweak one worsens the Higgs hierarchy problem [7,8].
Nevertheless, the extreme lightness of neutrino masses may also be naturally explained through a symmetry argument [9][10][11][12]. Indeed, the Weinberg operator violates B − L, which is an accidental and non-anomalous symmetry of the SM. This fact is exploited by different low-scale variants of the Seesaw mechanism, such as the inverse [13,14] or linear [15,16] Seesaws, to lower the mass scale of right-handed neutrinos. This way, the Higgs hierarchy problem is not worsened, and these new "Heavy Neutral Leptons" (HNLs) may be probed at laboratory energies. Still, even if light, the interactions of the HNLs with the rest of the SM particles are also suppressed by an unknown but small mixing with active neutrinos. The HNLs then constitute long-lived particles with "weaker-than-weak" interactions.
In this context, if additional new physics is present above the electroweak scale, it can have an important effect not only on SM processes but also on the behaviour of HNLs, possibly dominating over its weaker-than-weak interactions. To explore the impact of these high-energy new physics scenarios on HNLs, we turn to the framework of effective field theories (EFTs).
Effective field theories are an extremely powerful and versatile tool that encodes the indirect effects of new physics appearing at energy scales above a given physical process. They rely on an expansion in inverse powers of the new physics scale, introducing operators of dimension higher than 4; the higher the dimension, the more suppressed the operator is. In this context, the naive expectation is to see the first evidence of new physics from the Wilson coefficients corresponding to the lowest dimension and least suppressed operators. Thus, it is very suggestive that the only d = 5 operator that can be built out of the SM particle content is the Weinberg operator, which can account for the evidence of neutrino masses and mixings. EFTs are renormalizable order by order, with a finite number of counterterms required at each order of the EFT expansion. Therefore, all possible operators that can be built with low-energy particle content which are compatible with the symmetries of the theory should be included at each order. In this context, the Standard Model Effective Field Theory (SMEFT), that is, the series of higher dimension operators that can be built from the SM particle content, is a very useful and successful tool to constrain generic new physics (see, for instance, Ref. [17] for a recent review).
Since the evidence for non-zero neutrino masses can be explained by HNLs light enough to be included in the low energy spectrum, it is then interesting to extend the SMEFT description with a basis of new operators also incorporating the HNLs as fundamental building blocks: the νSMEFT .
In this work, we explore the effective operators (up to d = 6) that form a basis including these new particles and study the leading constraints on them from different sets of existing observables. For the first time, we present a comprehensive list of bounds that apply to each operator as a function of the HNL mass. As a first approach, we will only consider one operator at a time. We leave a more realistic global analysis for future work. Nevertheless, this approach generally leads to conservative constraints, as it does not allow for HNL production and detection through different operators. When relevant, we will comment on exceptions to this rule when cancellations among different operators lead to flat directions. In the same spirit, we will also assume negligible mixing between the HNLs and the active neutrinos. Mixing can lead to new production mechanisms or decay channels that would strengthen the constraints.
This paper is organized as follows. In Section 2, we introduce some generalities about HNLs, summarizing the main constraints on their mixing in Section 2.1. Then, in Section 3, we briefly introduce the νSMEFT and the new operators that extend the SMEFT by including HNLs. In Sections 4 to 7, we describe the different types of dimension-6 (d = 6) operators, the observable effects they can produce, and how current limits can constrain their Wilson coefficients. We present our conclusions in Section 8.

Heavy neutral leptons at the renormalizable level
We will extend the SM particle content by only one 1 right-handed neutrino field N ′ , singlet under the SM gauge symmetry group. The renormalizable part of the Lagrangian reads where L stands for a left-handed lepton doublet of flavour α, M N is the Majorana mass, and y α are the Yukawa couplings of the HNL to the different flavours of lepton doublets. The extension of the neutrino sector implies that the flavour eigenstates will generally contain some fraction of the heavy component due to the mixing U αN = y α v/ √ 2M N : where ν α(i) denote the interaction (mass) eigenstates and N is the heavy, almost sterile, neutrino. If the mixing elements U αN → 0 for all α, then N can be identified with the flavour state N ′ and M N is the physical mass of the heavy state. In this work, we assume the mixing elements are small and organize our discussion around the physical field N . As Eq. (2.2) shows, HNLs would participate in any process in which active neutrinos appear, "inheriting" their interactions with an additional suppression from the mixing. These elements are typically small, giving rise to weaker-than-weak interactions, long lifetimes, and suppressed production. Nevertheless, laboratory experiments can be sensitive to HNLs in a wide mass range.
Experimental searches have not found any significant evidence for their existence, but have provided strong limits on the mixing elements, U αN . These bounds are flavourdependent, as the processes relevant to HNL searches are considerably different depending on the flavour structure of the mixing matrix [55][56][57][58]. Constraints are typically stronger on the mixing with the electron and muon flavours, and less so on the mixing with tau. The standard assumption followed by experimental collaborations when providing limits on either U eN , U µN , or U τ N (strictly speaking, on their moduli squared, which are the only observable quantities) is to consider a simplified flavour structure. It is assumed that the mixing with one flavour dominates over the other two elements, whose impact on the phenomenology is assumed to be negligibly small. In this work, we will follow this simplified scenario. However, these commonly used benchmarks do not reflect the flavour structure typically found in realistic neutrino mass models. As a better approximation to the phenomenology of neutrino mass models, two additional benchmarks have been recently proposed in the context of the CERN's Physics Beyond Colliders initiative [59].
A comprehensive collection of experimental limits on HNLs is available on github . This database adds to the existing collection of limits on other light particles, such as the limits on axions and dark photons [60,61], as well as complementary efforts for HNLs [62].

Constraints on heavy neutrino mixing
Many bounds on effective operators, studied in the following sections, will be obtained by reinterpreting those existing on HNL mixing. Thus, we will first summarize the most relevant and up-to-date constraints existing in the literature. For previous studies summarizing constraints on HNLs, see Refs. [62][63][64][65][66][67][68][69] as well as recent reviews [70,71]. As mentioned, single flavour dominance is assumed, so we will discuss limits that apply separately on either |U eN |, |U µN |, or |U τ N |. Many experimental bounds are available, but we will only mention the leading limits in each HNL mass window here.
Electron flavour dominance (|U eN |): At the lowest masses, the best limits are provided by searches for an invisible resonance in π + → e + N decays. The best constraints are given by these "peak searches" at TRIUMF [72] and PIENU [73], including phenomenological uses of their data [74]. Below 16 MeV, Borexino has searched for the N → νe + e − decays of HNLs produced in 8 B decays in the Sun, providing competitive limits [75]. Above the pion mass, the dominant constraints come from peak searches in leptonic kaon decays, K + → e + N , in particular, those performed by the NA62 collaboration [76]. HNL decayin-flight signatures at neutrino experiments in a similar mass region provide competitive constraints with a complementary technique. The T2K collaboration searched for HNL charged-current (CC) decay channels in their near detector, ND280, using HNLs produced in kaon decays at the target [77]. Above the kaon mass, the CHARM [78] and BEBC [79,80] experiments set the strongest limits. Finally, for HNLs with masses above ∼ 2 GeV, HNLs are no longer copiously produced by meson decays at beam dumps, and colliders provide the best environment for detection. One of the most competitive limits was provided by the DELPHI collaboration [81] studying Z decays to produce the HNLs, and searching for their prompt or displaced decay inside the detector. Because the production proceeds through the neutral-current (NC) coupling, this bound is largely insensitive to the flavour structure of the HNL mixing. At the LHC, the ATLAS [82,83] and CMS [84] collaborations have set some of the strongest limits by using both W or Z boson decays and searching for the subsequent HNL decays into dilepton final states. Finally, for masses beyond the reach of collider energies, the HNL mixing can only be probed indirectly by testing the unitarity of the 3×3 PMNS matrix, which involves only the 3 light mass eigenstates and active neutrino flavours, through flavour and electroweak precision data [68,85].
Muon flavour dominance (|U µN |): This case follows a very similar pattern to the one above. As the muon mass is close to that of the pion, HNLs can only be produced in pion decays if they are lighter than about 30 MeV. Just as before, peak searches at PSI [86] and the PIENU experiment [87] provide the best limits. Above the pion mass, the decayin-flight searches at the PS191 and T2K experiments provide the strongest bounds [77,88,89] 2 . The region in between m π − m µ ∼ 34 MeV and m π was not targeted by any dedicated experimental search, but phenomenological studies have recasted limits from MicroBooNE [92] and T2K [90], as well as from measured kaon decays (KEK) [93,94]. Above M N ∼ 180 MeV, the E949 and NA62 kaon factories provide the strongest constraints using peak searches in leptonic kaon decays [95,96]. Above the kaon mass, low-energy neutrino experiments and kaon factories become insensitive, and high-energy beam dumps dominate by looking for the decays of HNLs produced in D meson decays and muon-neutrino upscattering in the detector. The leading limits are provided by BEBC [97], NuTeV [98], and CHARM [78,79]. Finally, for HNLs above ∼ 2 GeV, high-energy colliders stand out. The ATLAS [82,83] and CMS [84] collaborations provide the most stringent constraints using displaced decays of HNLs. Finally, DELPHI [81] also provides competitive limits using prompt or displaced decays from HNLs produced in Z decays. As for the electron case, for the highest HNL masses constraints can be set through flavour and electroweak precision data probing the leptonic mixing matrix unitarity [68,85].
Tau flavour dominance (|U τ N |): Finally, the mixing with the tau flavour is the less probed. Solar neutrinos can produce light HNLs (M N < 16 MeV) inside the Earth, and their subsequent decays in large-volume neutrino detectors, such as Borexino, provide the best limits [99]. A similar strategy was pursued by Ref. [100] to obtain bounds from Super-Kamiokande on HNLs produced by atmospheric neutrinos scattering in the Earth. Atmospheric neutrino data [101] provide a limit of |U τ N | 2 < 0.13 (90% C.L.) [102] for HNL masses in the 10 eV ≲ M N ≲ 10 MeV range [103]. For larger HNL masses, only high-energy beam dumps can produce them, through D and τ decays, so the best limits are those obtained in phenomenological reanalyses of CHARM [104,105] and BEBC [79,80] data. A dedicated experimental search was also performed by the ArgoNeuT collaboration [106]. While peak searches in π and K decays are not possible in this case, a peak search in τ → N 3π decays was performed by the BaBar e + e − collider [107]. As usual, for HNLs above the GeV scale (roughly 1.5 GeV), the most competitive limits are given by colliders. As explained above, the production of HNLs at DELPHI was mostly insensitive to the flavour structure of mixing, so DELPHI provides the strongest constraints [81]. As for the other flavours, the constraints for the highest HNL masses stem from unitarity bounds [68,85]. Figure 1 displays all the bounds discussed above, as a function of the heavy neutrino mass, for each lepton flavour dominance case. The constraints set by each experiment are displayed independently. It is worth mentioning that, besides these laboratory constraints, HNLs mixing with active neutrinos would have a considerable impact on the cosmological evolution of the Universe. Indeed, even though their mixing may not be large enough to thermalize HNLs with the SM plasma, they would still be produced, and, upon becoming non-relativistic, their contribution to the energy density of the universe could be sizable. For small masses and mixings, the lifetime of the HNLs could be larger than the age of the Universe, and their predicted population could exceed the measured dark matter component. Conversely, for heavier masses, the HNL decays could lead to an unacceptably large contribution to the energy density in radiation, N eff . Their decay products to charged particles could also raise the ionization floor after recombination and alter the CMB temperature and polarization power spectra. The combination of these effects allows setting very stringent bounds, stronger than the laboratory constraints summarized in Fig. 1 in some parts of the parameter space, as shown in Refs. [108,109]. Similarly, if HNLs decay during BBN, they might alter its predictions for the primordial abundance of light elements, leading to stringent constraints [110][111][112][113][114]. For larger mixing angles, however, HNLs decay before the BBN epoch, and cosmological constraints are no longer effective. In that sense, they provide an upper as well as a lower limit for the mixing of HNLs, with a direct complementarity to laboratory-based bounds.
The aforementioned constraints also rely on a given cosmological model and its extrapolation to the early universe, and it is fair to treat them separately compared to the more direct laboratory limits. The two sets of bounds should therefore be regarded as complementary since a disagreement between them could signal a deviation from the standard cosmological model. Thus, in this work, we will focus only on the more direct laboratory and astrophysical constraints, but the existence of these complementary cosmological bounds should also be acknowledged.

Including HNLs in SMEFT
The Standard Model Effective Field Theory (SMEFT) is a very powerful tool to parametrize the effects on low-energy observables of new physics lying at high energies. The possible BSM effects are encoded in effective, non-renormalizable operators of dimension 5 and higher, which are built from the SM fields and respect its fundamental symmetries. Thus, the SMEFT Lagrangian consists of an infinite tower of operators: where Λ is the scale of new physics, up to which the effective theory is applicable, and O is the dimension d = N effective operator. Here, C i are the Wilson coefficients of each operator: dimensionless parameters that control their respective coupling strength. The higher the dimension, the greater the number of possible operators, but the smaller their effects on low-energy observables, as they are increasingly suppressed by larger powers of the new physics scale. The effective field theory framework is not suitable if the new particles are lighter than the energy scales of the observables. In that case, the new propagating degrees of freedom have to be added to the theory. This is the case for the light HNLs that we consider since M N can be much smaller than the electroweak scale. With the addition of the new field N , a singlet under the SM symmetry group, the SMEFT framework needs to be extended. This new framework will contain new operators involving the field N and other SM particles. This is a compelling approach to thoroughly explore the potential effects of new light particles and takes advantage of the power of EFTs to describe generic and indirect new physics effects.
The operators of this extended theory, which we refer to as νSMEFT, were first systematically considered in Ref. [20] (see also [18,19]), and then in Ref. [27], where redundant operators were removed, and new dim-7 operators were included. Many of these nonrenormalizable operators provide a new way for HNLs to interact with other SM particles. In addition, many processes mediated by the effective operators are also generated by the renormalizable interactions described in Section 2.1, namely, the mixing between active and heavy neutrinos. Thus, the bounds arising from this kind of observable may also provide limits on the corresponding Wilson coefficients. Whether the existing limits apply and in which way depends on the particular operator considered.
Our main focus will be on d = 6 operators. At higher dimensions, the operators are further suppressed and their number quickly gets out of hand, hindering the extraction of conclusions from an EFT. Let us briefly comment on the only two operators that exist at dimension 5 (see Ref [29] for a detailed discussion), 3) The first one contributes to the Majorana mass term and to the invisible decays of the Higgs boson into two HNLs. The associated phenomenology was studied in [19,29,33,37,49]. The focus of these works was on the sensitivity of collider experiments to prompt or displaced vertex signals originated by HNL visible decays, which depend on the value of the mixing U αN and the HNL mass. Additionally, the ATLAS [115] and CMS [116] measurements of the Higgs boson signal strength, µ, provide an important constraint on C d=5 Higgs /Λ, since the exotic Higgs decay H → N N would affect this quantity. Following [117] and using the combined bound µ ≥ 0.94 (95% C.L.), we find C d=5 This constraint is comparable with the expectation for the LHC sensitivity from direct searches [29] or even stronger. It should be noted that this current bound is independent Decay-in-flight and peak searches for e and µ. PMNS unitarity and peak searches for τ . Moments O α NB L α σ µν N HB µν Neutrino upscattering, monophoton searches.
Monolepton searches, decay-in-flight and peak searches. Table 1. The dimension-six operators of the νSMEFT involving at least one N field. We list the leading constraints on each one of the operators in the range 1 MeV < M N < 100 GeV. We also include the bound on the dimension-five operator involving HNLs and the Higgs derived in this work.
of the HNL mixing and mass as long as M N < 40 GeV. The bound becomes increasingly weaker for heavier masses up to M H /2. The second operator requires at least two different HNL fields to not vanish and amounts to a transition magnetic moment for the HNLs involved. It has been discussed, for instance, in Refs. [21,47].
In the following sections, we present all the d = 6 operators involving at least one N field. We will discuss them one at a time, summarizing the processes they mediate and how they are constrained to finally display the bounds on the corresponding Wilson coefficients as a function of the HNL mass. Barring flat directions, which we will try to comment on, this approach is generally conservative and provides the first step towards a more general and ambitious global fit in the future. Details on the procedures followed to derive and recast the constraints are given in Appendices A and B. A list of the operators and the corresponding limits is shown in Table 1.

Higgs-dressed mixing
The first and simplest operator we consider is the addition of a Higgs bilinear to the Yukawa coupling of Eq. (2.1), This operator is a neutrino Yukawa coupling with an extra pair of Higgs doublets. Once they acquire vevs, this operator induces extra contributions at tree level to light neutrino masses and their mixing with the HNLs. In the limit where M N ≪ v ≪ Λ, the generated mass and mixing read Note that there are three copies of this operator, one for each lepton flavour, whose Wilson coefficients are, in principle, independent. At low energies, the mixing induced by this operator gives rise to an identical phenomenology to that already studied in the case of the Yukawa coupling of Eq. (2.1); thus, all bounds on mixings from direct constraints apply directly, yielding limits on the Wilson coefficient. There could be possible cancellations between this contribution and the tree level d = 4 Yukawa couplings, leading to possible flat directions; as we work under the assumption of negligible standard mixing, we will not take them into account. In principle, it could be possible to lift this flat direction, as this effective operator also mediates processes involving two HNLs and several Higgs bosons, and the d = 4 Yukawa does not. However, these processes, such as HH(H) → N N , require multi-Higgs production and thus do not offer promising experimental prospects.
The existing limits on the absolute neutrino mass scale would also allow setting extremely stringent constraints on the Wilson coefficient. However, the expression shown in Eq. (4.2) corresponds to the inclusion of a single HNL. When several of these new particles are considered, their contributions to neutrino masses could cancel each other, while this is not a possibility for their mixings. Indeed, this is exactly what would be expected in low-scale realizations of the Seesaw mechanism. In those scenarios, the contribution to the neutrino mass cancels, due to its protection from the approximate B − L symmetry, while the mixing does not. Such a protection mechanism needs to be in place to have such light HNLs in the first place. Thus, we do not include bounds on the Wilson coefficient arising from the current information on the absolute value of neutrino masses.
Apart from contributing to the neutrino masses and mixing, this operator induces an invisible Higgs decay channel, into a light and a heavy neutrino. The direct limits on the invisible branching ratio of the Higgs are B(h → N ν α ) < 0.13 , as reported by ATLAS [115], and B(h → N ν α ) < 0.16, as reported by CMS [116], both at 95% C.L. In our work, we will combine the most stringent upper limits based on the Higgs signal strength by ATLAS [115] and CMS [116], yielding B(h → N ν α ) < 0.06 at 95% C.L. following the Feldman-Cousins [118] prescription.
In Fig. 2, we show the bounds on the Wilson coefficients of this operator under the three lepton-flavour dominance scenarios, shown in three separate panels. Most shaded regions (except the dark grey bands) correspond to reinterpreted mixing bounds, and, according to the assumption of single flavour dominance, are different for each lepton flavour. The limits from Higgs invisible decay (dark grey region) apply to the dominant Wilson coefficient, independently of its flavour.

Bosonic currents
It is possible to couple bosonic currents to the V + A currents constructed out of heavy neutrinos. The resulting operators read where ℓ α is a right-handed charged lepton of flavour α. After electroweak symmetry breaking, the first operator will yield a coupling of the Z boson to two HNLs, while the second operator will produce a coupling of the W boson to one HNL and a charged lepton. The latter will exhibit three copies, depending on the flavour of the charged lepton. The phenomenological consequences of these two operators are quite distinct, and so are the constraints placed upon their corresponding Wilson coefficients. They are discussed separately below.

Neutral bosonic currents
The operator O HN yields a vertex between two HNLs and a Z, namely C HN Λ 2 gv 2 2c W N γ µ N Z µ , opening a new decay channel for this boson if kinematically allowed by the masses of the HNLs. This decay, controlled by the Wilson coefficient C HN , constitutes an extra contribution to the invisible decay of the Z boson, which was measured with high accuracy at LEP [119,120]. We have performed a χ 2 fit to this observable, obtaining upper limits on C HN /Λ 2 .
Similarly, this neutral bosonic current would mediate decays of neutral mesons (π 0 , η, η ′ , J/ψ, ρ 0 , Υ(1S)) into two HNLs, if light enough. Analogous invisible decays are also present in the SM, with light neutrinos in the final state [121,122]. For pseudoscalar mesons, these processes require a chirality flip, so in the SM they are suppressed by the small neutrino masses; for vectors, they are very subdominant with respect to hadronic and electromagnetic channels. Thus, this kind of decay has not been observed so far, and only upper bounds are available. Imposing the decays into HNLs to respect the experimental limits provides another method to constrain the Wilson coefficient. We will only display the upper bound on the neutral pion and the Υ(1S) invisible decay, as they are the most stringent, depending on the HNL mass. In particular, the NA62 collaboration determined the corresponding branching ratio of the pion to be smaller than 4.4 × 10 −9 , at 90% C.L. [123]. For the Υ(1S), we employ the limit of 3 × 10 −4 , at 90% C.L. [124]. Stronger limits on this invisible branching ratio can be obtained using γγ → π 0 → N N in supernovae (see also the supernova discussion below). The limit is 3.2 × 10 −13 for the BR [125].
Another probe of this effective operator comes from collider physics, in particular, monophoton searches performed at LEP [126]. This channel, intended to look for processes in which one or more invisible particles were produced in the e + e − collision, has the photon as the only detectable signal. These searches are blind to the nature of the invisible particle, dark matter being one of the most recurrent candidates. We follow Ref. [127], where an EFT framework was also adopted, with dark matter fermion singlets playing a similar role as our HNLs. We translate this information into constraints on C HN /Λ 2 . Note that Ref. [127] employs a basis of effective operators different from the one considered here. To derive our bounds, we compute the monophoton production cross section for the operators of our basis, applying the same cuts in the photon energy and angle, and rescale the constraints accordingly to obtain the same cross section in the signal region as for the operators in Ref. [127]. Details on this rescaling procedure can be found in Appendix A.3.
A very similar process could be mediated by this operator replacing electrons with quarks and the photon with a gluon, yielding experimental signals composed of a single jet and missing energy. Thus, monojet searches at the LHC pose an additional constraint on the corresponding Wilson coefficient. Limits from these processes are not available for the operator under consideration, and a proper rescaling would involve a detailed simulation, which is out of the scope of this work. Nevertheless, a crude estimate, obtained by rescaling parton-level cross sections, yields bounds similar to those from monophoton searches. These searches will also be relevant for the observables discussed in Section 7.
Supernovae also offer competitive probes of particles that couple very feebly with the SM. Due to their weakly interacting nature, neutrinos are the dominant cooling mechanism of core-collapse supernovae. This mechanism is compatible with the observation of neutrinos from SN1987A by the Kamiokande-II [128] and IMB [129] neutrino detectors. Thus, a BSM particle would be incompatible with observations if it exhibited a mass and couplings such that it extracted energy from the supernova at a faster rate than SM neutrinos [130]. This usually provides upper and lower limits on the couplings: if the new particle couples too strongly with the SM, it would be trapped inside the supernova, whereas, if the coupling is too small, the new particle is produced much less often than neutrinos, cooling the supernova less. We follow Ref. [131], where these arguments are used to constrain the parameter space of a fermion singlet, also working in an EFT framework. We again rescale their results so that the production and scattering cross sections of the HNLs coincide with those for the operators in Ref. [131] along the lines defining the lower and upper bounds respectively (see Appendix A.2). Note that, depending on the HNL mass, the most stringent upper bound may be given by HNL production in π 0 decay in supernovae, as discussed above.
Finally, it should be noted that for this operator, or, more generally, for those involving two HNL fields, and in absence of mixing or other interactions, HNLs would be completely stable. Thus, their production in the early Universe could lead to an overabundant dark matter component and very stringent constraints. As discussed at the end of Section 2, we regard these important cosmological constraints as complementary to the direct bounds discussed in this work. Indeed, their comparison provides a test of the necessary assump- tions on the cosmological model, which are not present in the laboratory or astrophysical observations we consider here. For instance, if the reheating temperature were below the HNL mass, the overproduction constraint would not apply. Figure 3 collects the bounds on C HN /Λ 2 derived from the different observables discussed above: invisible π 0 , Υ(1S), and Z decays (white, yellow and dark blue regions respectively), monophoton searches (pink area) and supernova cooling arguments (light blue region). Note that the bounds on the standard HNL mixing do not apply, as they typically involve CC interactions (mostly charged mesons decay), which are not mediated by this operator.

Charged bosonic current
After EW symmetry breaking, the operator O α HNℓ leads to a vertex between an HNL, a charged lepton, and a W , with a similar chiral structure to the SM charged currents: Thus, this operator yields a coupling to the W to a right-handed current, analogous to the left-handed counterpart inherited through mixing, controlled by the Wilson coefficient: We then reinterpret the current bounds on HNL mixing as limits on C α HNℓ /Λ 2 . Most of the bounds only involve processes mediated by the SM charged currents, and thus translate directly by means of U CC αN = U αN . However, some of the experimental searches also consider NC-mediated processes, absent for the effective operator under consideration. The limits provided by those searches need to be rescaled, to remove the neutral contributions while keeping constant the expected number of events at the detector. We describe our rescaling procedure in Appendix A.1.
Supernova cooling and LEP monophoton searches are also able to probe charged bosonic currents. There are no basic conceptual differences with the neutral case; the main distinctions appear at a diagrammatic level, as the HNLs are now produced via the exchange of a W in a t channel. However, the associated bounds are not competitive with those arising from the mixing, as the corresponding processes are NC-like. Thus, we will not display them.
Finally, this effective operator mediates leptonic decays of muon and tau, ℓ α → ℓ β ν β N . We have performed a χ 2 fit of the corresponding decay widths of the muon and the tau to the experimental measurements [119], finding limits on the Wilson coefficient that dominate at low masses. Note that, in the channel τ → µνN , the experimental determination exhibits a slight tension, of roughly 2σ, with the SM prediction. This means that, at this confidence level, the Wilson coefficient in the tau flavour dominance scenario is constrained to a band instead of only an upper bound. We display the corresponding lines in a dashed style, as they are not completely equivalent to those arising from other observables. Figure 4 contains the bounds on the Wilson coefficients of this operator, coming from reinterpretations of limits on heavy neutrino mixing. Once again, we assume single flavour dominance, displaying on different panels the relevant constraints for each lepton flavour.

Tensor current
Similar to the d = 5 operator in Eq. (3.3), dipole couplings between the HNL and the lepton doublet can be considered. Two different operators arise: Three copies of these two operators exist, one for each lepton flavour. Once again, their Wilson coefficients are, in principle, independent. After the Higgs develops a vacuum expectation value, these operators are translated into three dipoles for the HNL, coupling it to the photon, the W , and the Z: Bounds on these dipoles can be translated into limits on the Wilson coefficients by applying the relations above. Most limits on the Wilson coefficients C NB and C NW are obtained from processes mediated by the magnetic dipole moment d γ . The weak dipoles d Z and d W are less important for low-energy processes but can contribute to HNL production at colliders, for example. We recast the LHC and LEP limits obtained in Ref. [132] to match our one-operator-at-a-time approach 3 . We note that a clear flat direction is present for C NB /C NW = − tan θ W so that the dipole moment with the photon would vanish even with large Wilson coefficients. In this case, none of the limits shown in Figs. 5 and 6 apply, as they rely either on N production or decay via d γ . We leave a study of this particular case to future literature. [133][134][135][136] or HNL production and decay [99,134]. The LEP monophoton searches and supernova cooling bounds, mentioned in the previous sections, also constrain the dipole moment. Finally, this quantity must also be small enough to respect primordial abundances after BBN [133]. See [137] for a recent review on all these bounds. However, as we have already discussed in Section 2.1, these relevant bounds ultimately depend on the cosmological model, and we regard them as complementary to the laboratory constraints studied here.

Probes of magnetic dipole moments include searches for exotic electromagnetic interactions of light neutrinos
In the scenario in which the couplings lie along the flat direction, couplings to the W and/or Z must be considered instead. In principle, the operator O NW would also mediate CC interactions, such as the ones involved in the standard bounds on the mixing U αN , but with a different Lorentz structure. Similarly, O NB would induce couplings between HNLs and neutral mesons. Nevertheless, the couplings with mesons vanish due to their particular Lorentz structure. Thus, all bounds arising from HNL production or decay via meson interactions do not apply to these operators either.
Finally, the operator O NB would mediate a new invisible Z decay channel, into an HNL and a light neutrino. Imposing its rate to be consistent with current measurements on invisible Z decay would yield yet another limit on C NB /Λ 2 . However, we find this bound to be an order of magnitude worse than the one arising from monophoton searches, and do not include it in our plots. Figures 5 and 6 show the present bounds at 90% C.L. for C NB /Λ 2 and C NW /Λ 2 , respectively. Both of them contain three panels, one for each lepton flavour, following the single flavour dominance assumption. Most bounds on the dipole moments depend on the flavour of the involved neutrino, so the different panels will, in general, contain limits provided by different experiments. Some sources of constraints are flavour independent, such as those provided by LEP or by supernova cooling arguments.

Four-fermion interactions
Several different four-fermion operators involving HNLs can be written at d = 6. We will separate their study into two main sections, depending on whether the operators mediate CC-like or NC-like processes.
Analyzing their potentially complex flavour structure leads to a large number of operators of each type, a priori independent. As in the previous sections, we will consider separately the operators that affect electrons, muons, and taus, treating the corresponding Wilson coefficients as independent parameters, to be bounded by different observables.
To derive the constraints on the effective operators involving quark fields, their flavour structure becomes very relevant. Indeed, one of the most advantageous ways of producing HNLs is to exploit decays with a relatively low rate in the SM to compete against, such as K or B meson decays, which must proceed through CKM-suppressed weak processes. Thus, these searches only apply if quark flavour violation is allowed in the effective operators.
However, some operators contain both charged and neutral couplings and, if allowed to be flavour-violating, could lead to much more strongly constrained flavour-changing neutral processes, such as neutral meson oscillations at the loop level. In the spirit of deriving conservative limits on the Wilson coefficients, we will assume that the effective operators are flavour universal and diagonal in the SM flavour basis. Thus, the only source of flavour violation in the quark sector will be the SM CKM matrix, which will appear in all CC-like operators, in an analogous way to the SM CC interactions when going to the mass basis. Indeed, we will also assume that the corresponding rotation for the right-handed fields, which will be physical for some operators, also corresponds to that of the CKM matrix 4 .
For comparison, for the flavour-changing CC-like operators, we will show our results both under the assumption of a CKM suppression and for the flavour-blind case, where all flavours couple with the same intensity (all flavour copies share a common Wilson coefficient). In this way, it will be apparent which parts of the bounds derive purely from the experimental results, and which rather stem from the flavour alignment prescription. However, notice that, in the flavour-blind case, stronger constraints from FCNC processes could generally apply; we only display this case for the sake of comparison.

Neutral currents
Three sets of four-fermion operators mediate NC-like processes. The first one reads where f can stand for any right-handed fermion in the SM, either leptons or quarks. In the latter case, we will consider couplings to the up-and down-type quarks independently, through the coefficients C uu and C dd . For simplicity, we assume that this coefficient is generation independent. Analogously, very similar operators can be written out of left-handed SM fields, where Q i denotes a quark doublet of flavour i. In the case of leptons, we will once again treat the three flavour copies of O LN as operators with independent coefficients. In the case of the quarks, we assume that the coefficients are generation independent. Finally, it is also possible to construct an independent operator with a scalar structure: The fact that none of these operators mediates CC-like processes makes them difficult to probe. For instance, most bounds on HNL mixing cannot be reinterpreted for this purpose, as they rely on the production of HNLs in the decays of charged mesons. Even the constraints from DELPHI, based on HNL production in Z decays, cannot readily be used. Nevertheless, O uu , O dd and O QN contribute to invisible decay channels for neutral mesons, with the π 0 and the Υ(1S) being the most stringent. These operators can also induce quark scatterings in which the only visible signal (aside from the HNLs) is a single jet, produced mainly by a gluon emitted by any of the quarks. Such monojet searches have been performed in the LHC. Ref. [31] recasted the experimental limits into bounds on the corresponding Wilson coefficients. As assumed in [31], the HNL mass can be neglected in this process for the range of masses we explore, corresponding to a straight line in our plots.
We find no observables able to constrain the operators involving two muons or two taus. In the case of two electrons, LEP monophoton searches and supernova cooling arguments apply analogously to previous sections, as electrons are present in the possible production of HNLs. As discussed before, a rescaling is required to account for the different cross section each operator would lead to (for details, see Appendix A.2). The operators O ee and O LN have different chiral structures; however, the bounds on their corresponding Wilson coefficients turn out to be the same. The same observables also constrain C αβ LNLℓ , provided α = β = e. If α ̸ = β, ℓ α → ℓ β νN decays are mediated by this operator. We perform a χ 2 fit of the corresponding rates to the measured µ and τ leptonic decays, obtaining bounds for C eµ LNLℓ , C eτ LNLℓ and C µτ LNLℓ . Figures 7 to 9 show bounds, at 90% C.L., on the mentioned Wilson coefficients. The main limits arise from invisible meson decay and monojet searches in the case of quarks, monophoton searches and supernova cooling for electrons and muon and tau decays for the other lepton flavours..

Charged currents
Three types of four-fermion operators appear at d = 6, mediating CC-like interactions. The first one involves a right-handed charged lepton, a right-handed down-type quark of flavour i, and a right-handed up-type quark of flavour j, namely where u i and d i are the right-handed quark fields and i is the generation index. In the flavour-blind case, Z duNℓ ij = 1, while, under the flavour alignment hypothesis, an insertion of the corresponding CKM matrix element V ij would be needed: It is also possible to write a scalar coupling of the form As before, under the flavour alignment hypothesis, we assume this operator is flavour universal and diagonal in the flavour basis, and, upon rotating to the mass basis, the CKM mixing matrix will control the degree of flavour violation. while in the flavour-blind scenario Z LNQd ij = 1. Note that exchanging the down-type quark and the HNL fields leads to a different operator, which shares the flavour coefficients of the previous operator. However, through a Fierz transformation, it is possible to rewrite this operator as so that the bounds on the previous operator will also apply to this one. Finally, the last independent operator takes the form (7.12) in the case of flavour alignment, and Z QuNL ij = 1 in the flavour-blind scenario. The set of operators in Eqs. (7.5), (7.7) and (7.11) mediate charged currents involving quarks and HNLs, producing interactions between charged mesons and HNLs. These processes are the main source of bounds on the standard HNL mixing, so the corresponding limits can be translated into constraints on the Wilson coefficients. The interactions mediated by the effective operators exhibit different Lorentz structures than the SM, so they yield HNL production and decay rates different from those mediated by standard mixing. Thus, as for the previous analyses, a rescaling of the existing bounds is necessary, to account both for the absence of neutral currents and for the different production and decay rates of the HNL. This procedure was recently advocated and applied to this type of operator in Ref. [52] for some example observables.
All these operators could also induce monolepton processes in colliders, in which the final state consists of a single observed lepton and missing energy (carried by the HNL). Searches for these signatures have been performed at the LHC and were recasted to constrain the Wilson coefficients of these operators in Ref. [31]. We directly employ their bounds, which are independent on the mass of the HNL (as the typical energies at the LHC are much larger than the range of masses under discussion).
We will not discuss the phenomenology associated with O QuNL , as it is identical to that of O LNQd . This can be inferred from Eq. (7.10). The term with a tensor structure does not contribute to interactions with mesons. On the other hand, the term with a scalar Lorentz structure is identical to O LNQd (modulo a factor of 2), thus yielding no different phenomenology. The bounds on C LdQN can be directly read from those on O LNQd by rescaling them by a factor of two.
The operators O LNQd , O LdQN and O QuNL also mediate NC-like processes involving an HNL, a light neutrino, and a qq pair. This would induce invisible decays of neutral mesons. Imposing the experimental limits on those processes could also constrain the corresponding Wilson coefficients. However, these bounds are always considerably looser than those arising from charged currents. As the limits would apply to the same Wilson coefficients, we will focus on those derived from CC-like processes.   by the dotted lines with lighter colours, flavour blindness is assumed (Z ij = 1), whereas, in the regions delimited by the solid lines with darker colours, we adopt the flavour alignment prescription.

Conclusions
The addition of right-handed neutrinos, or HNLs, to the SM particle content, is arguably the simplest way to explain the evidence for neutrino masses and mixing. The mass scale of these new particles is thus a new dimensionful parameter of the model, to be determined empirically. If this mass scale is below EWSB, the HNLs should be considered, together with the other SM fields, as fundamental building blocks of any EFT aiming to probe new physics at a higher energy scale.
In this work, we have systematically discussed the operators that would extend the basis of the successful SMEFT paradigm with the extra HNL building blocks. This new EFT is referred to as νSMEFT. For each new operator in the νSMEFT basis, we have discussed the laboratory and astrophysical bounds that would apply in absence of any other new physics contribution. That is, we consider one operator at a time and neglect also any contribution from the mixing of the HNLs with the SM neutrinos. This is the first step towards a more thorough and complete exploration of the operator basis. Furthermore, this procedure is generally conservative, as it prevents the possibility of different operators reinforcing each other's contributions or having different new physics at the production and detection of a given process. Nevertheless, we have discussed a few instances where flat directions in the operator basis could avoid certain bounds.
It should also be noted that these constraints are complementary to the ones that can be derived from the impact of the HNLs in the cosmological evolution of the Universe. Indeed, interactions between the HNLs and the SM fields would lead to their production in the early Universe. Depending on their lifetime, this could lead to too large contributions to the matter abundance, or their decay products could violate the present bounds on the radiation energy density present during CMB. They could also affect the CMB and BBN predictions. Since, unlike the laboratory constraints discussed in this work, these constraints rely on assumptions on the cosmological model, they should be explored together as a consistency check of those assumptions.
We find that the same experiments that constrain the mixing of HNLs with the SM neutrinos, such as peak searches and beam dump experiments, as well as collider searches at higher HNL masses, can be easily reinterpreted, as suggested in Ref. [52], to provide the most stringent direct bounds on the Wilson coefficients of CC-like operators with only one HNL field. Conversely, operators with two HNL fields of NC type are significantly more difficult to constrain in absence of mixing with the SM neutrinos. We find that these operators can instead be bounded from constraints on invisible decay of the Z boson or neutral mesons, supernova cooling limits, and monophoton and monojet searches at colliders.
We have presented plots summarizing the most relevant bounds derived on each Wilson coefficient as a function of the HNL mass over a very wide range of values, spanning from the MeV to the electroweak scale, 5 and made them available on github . Table 1 summarizes the operators and the corresponding limits. These results provide an overview of the present constraints applying to each operator in the EFT basis, enlarged upon the addition of HNLs, for future searches and analyses to improve upon. They also provide a first step towards more ambitious and realistic analyses allowing the presence of several operators, as well as mixing between the HNLs and the SM neutrinos simultaneously. under Grant CEX2020-001007-S, during the Extended Workshop "Neutrino Theories", where this work developed. Some effective operators considered here induce interactions akin to those generated by mixing between heavy and active neutrinos. While the standard mixing allows all the processes that could be possible in the SM, by simply replacing a light neutrino with a heavy one, the new operators may only generate a subset of those. We discuss how to reinterpret constraints on the mixing as constraints on the Wilson coefficients of the d = 6 operators. A similar rescaling procedure was also advocated in Ref. [52].
Charged bosonic currents. The operator in Eq. (5.2) induces an effective CC interaction. In analogy to the standard mixing, in Eq. (5.3) we defined the variable U CC αN ≡ C α HNℓ v 2 √ 2Λ 2 , which can be related to the standard mixing, U αN , in an experiment-dependent way.
The bounds set by peak search experiments apply equally to U αN and U CC αN , as the HNLs are always produced via charged currents, and their decay is not observed. In this case, the rescaling is trivial. Constraints obtained from HNL production in Z boson decay, such as the ones by DELPHI, would not apply.
For decay-in-flight searches, HNL production is identical to the standard mixing case, as it takes place primarily through charged meson decays. 6 The decays, however, are modified, and we take into account that NC channels are no longer available. By equating the number of events, the limits on the effective mixing of Eq. (5.3), can be obtained by, where Γ mixing (N → X) is the decay width of the HNL into the signal final state X (e.g., X = νe + e − , e + π − , . . . ) induced by the standard mixing scenario, and Γ CC (N → X) is the one induced by the charged bosonic current. The hat indicates that the decay width is deprived of its |U αN | 2 mixing factor or EFT coefficient C 2 /Λ 4 , such thatΓ mixing = Γ mixing /|U αN | 2 andΓ CC = Γ CC /|U CC αN | 2 . An analogous rescaling applies to collider constraints from ATLAS and CMS, where the HNL is produced in W ± decays, and the decay channels vary depending on whether NC was considered in the original search. In both cases, the Wilson coefficient will then be obtained by applying consecutively this ratio and the rescaling in Eq. (A.1).
Note that in the case of τ flavour dominance, decay-in-flight searches are insensitive to this operator, as the only possible decays of the HNL involve a τ lepton in the final state, which is either kinematically forbidden or not directly searched for.
Charged current four-fermion interactions. In a similar fashion to the charged bosonic current, the operators in Eqs. (7.5), (7.7) and (7.11) induce a subset of the interactions contained in the standard mixing case. As these operators involve quarks, the only processes they can mediate are mesons decay into neutrinos and vice versa. In that case, searches at CHARM, DELPHI, ATLAS, and CMS do not apply, as they tag HNL decays to fully leptonic final states. Furthermore, the rates of HNL production and decay can be modified due to the different Lorentz structures of the operators.
To illustrate the rescaling procedure, let us take the operator in Eq. (7.7), namely Here, Z LNQd ij stands for the flavour coefficient. If flavour alignment is assumed, Z LNQd ij = V ij . As usual, the necessary rescaling is achieved by equating the hypothetical number of events between the standard mixing case and the νSMEFT operator case. In the case of a peak search in meson decays, only the production decay rates are relevant. If the decay is, say, a pion into an electron and a neutrino, we can obtain a bound on the Wilson coefficient by means of where the decay width in the numerator is the one induced by mixing and, in the denominator, the one induced by O LNQd . Here,Γ LNQd = Γ LNQd × Λ 4 /|C e LNQd | 2 . This relation can be easily generalized for other meson decays and operators.
In the case of decay-in-flight searches, the decay rates of the HNL need to be taken into account. For instance, let us consider the same operator as above, O LNQd , in an electron flavour dominance scenario. In a search for both N → πℓ and N → ℓℓν decays in flight from HNLs produced by kaon decays at the target, only the former channel is sensitive to C e LNQd . The limit on the Wilson coefficient can then be extracted as The generalization to other lepton flavours and operators is straightforward and takes into account the different mixing assumptions and decay widths.

A.2 Supernova cooling
Supernova cooling arguments provide upper and lower bounds in the coupling of light particles. If the new states couple too feebly to the SM, their production is not enough to cool the supernova efficiently. On the other hand, if their interactions with the SM are stronger, they would be trapped inside the supernova, keeping the energy inside the system. In the intermediate regime between these two cases, the new particles can escape from the supernova, effectively cooling it.
Ref. [131] constrains the parameter space of a dark matter particle with a four-fermion vector coupling to the SM electromagnetic current, As HNLs play an identical role to DM fermions in the cooling of the supernova, their bounds can be readily translated as bounds on our operators, upon suitable rescaling, where C i and σ i are, respectively, the Wilson coefficient of a particular operator O i and the cross section of a process mediated by this operator. The rescaling factor is then obtained To compute the lower bounds, we employ the cross section for the scattering of the new fermion off the SM particles present in the supernova, whereas to compute the upper bound we employ the cross section for the production of the fermion. The particular processes will differ depending on the considered operator.
HNLs are mainly produced in e − e + annihilations. Thus, O QN cannot mediate this channel, and the operators involving charged leptons can only do so for electron flavour dominance (their muon and tau flavour copies cannot be constrained by supernova cooling arguments). The upper bound on the Wilson coefficient of the corresponding operator is obtained by simply rescaling the HNL production cross section: Note that an alternative efficient channel for HNL production is γγ → π 0 → N N . In fact, supernova cooling arguments provide a strong bound on the invisible π 0 decay. Imposing this limit yields upper limits for the Wilson coefficients of O uu , O dd and O HN . In the latter case, the mentioned constraints will compete with those given by e + e − annihilation, whereas they are the only source of upper bounds for O uu and O dd . Note that the quark content of the operator O QN is orthogonal to that of the π 0 , so this operator is unable to mediate its decay into HNLs, and thus cannot be constrained by supernova cooling arguments. The scattering processes that may affect the HNL will depend on the effective operator under consideration. The vectorial operator considered in Ref. [131] where the sum in X comprises all the possible HNL scatterings mediated by the particular operator O i . Of course, all these cross sections are energy dependent. The temperatures inside a supernova, and thus the energies at which these processes take place, vary considerably as a function of the distance with respect to the centre of the supernova. The abundances of electrons, positrons, protons and neutrons also exhibit radial dependences, which determine how likely their interactions are. In order to capture these effects, we compute averaged cross sections (denoted by σ in the equations above), integrated over the temperature profile of the supernova and convoluted with the abundance of the involved SM particles:  where T (r) is the temperature and n(r) is the number density of a given particle at a distance r from the centre of the supernova. σ is the energy-dependent cross section of each process (mediated by the operator under study), in which each particle is assigned the energy corresponding to the temperature. We extract the supernova temperature and abundance profiles from Ref. [131].

A.3 Monophoton searches
At LEP, searches for e + e − collisions with a single photon recoiling against invisible particles were used to derive limits on dark matter particles and HNLs. Ref. [127] places limits on DM production through effective four-fermion operators using these searches. From a collider perspective, HNLs and DM behave identically, as missing energy, so this information also allows constraining some of our effective operators, which are able to mediate e + e − → N N γ processes. Such is the case of the operators O ee , O e LN , O HN , O e HNℓ and O ee LNLℓ . To recast their limits on our Wilson coefficients, we account for the different Lorentz structures of the operators. In particular, one of the operators considered in Ref. [127] is a vectorial four-fermion coupling, similar to the one in Eq. (A.4). The bounds on our operators are then related to those on the vectorial one by where i runs through the operators mentioned above. To compute this ratio, we calculate the e + e − → γN N cross section, as described below. The 2 → 3 phase space is parameterized by the energy of one of the HNLs, E N , the energy of the photon, E γ , the relative azimuthal angle between them, φ N γ , and the polar angle of the photon, θ γ . The total cross section is given by where |M| 2 is the corresponding matrix element and E CM = 200 GeV is the center-of-mass energy [127]. The amplitude is determined by each particular operator and can be expressed in terms of the products of the four-momenta of the particles involved. It depends on the relative polar angle between the photon and the neutrino, θ γN , and the polar angle of the neutrino, θ N . The former is fixed once the energies of both particles are determined, while the latter is geometrically determined in terms of the angles in the CM, The energy of the photon can vary from 0 to E max γ = (E 2 CM −4M 2 N )/2E CM , where the latter is obtained when the HNLs travel in the same direction, opposite to that of the photon. Finally, the maximum (minimum) energy the HNL can carry depends on the energy of the photon and is obtained by solving cos θ γN = −1 (cos θ γN = 1), which gives We integrate over the physical kinematical range considering the cuts employed in the analysis by LEP, requiring the photon energy to be at least a 6% of that of the electron and its polar angle to be at least 45 • and no larger than 135 • . Note that the case of the operator O ee LNLℓ is slightly different, as it gives rise to a final state with a light neutrino and an HNL, instead of two HNLs. However, the recasting procedure is very similar, with some differences in the kinematics described above. The expressions can be easily generalized by assuming that one HNL is massless.

B Decay widths
In this appendix, we provide explicit expressions for the relevant decay widths of SM particles decaying into HNLs and vice versa, for each of the operators considered.  Table 2. Couplings controlling the interactions between neutral mesons and HNLs mediated by NC-like four-fermion operators. s 0(8) and c 0(8) stand for the sine and cosine of the mixing angles θ 0(8) .
Higgs-dressed d = 5 operator. As discussed in Section 3, this operator yields an invisible decay of the Higgs into two HNLs, given by the rate where M H is the Higgs mass, v its vev and x N ≡ M N /M H .
Higgs-dressed mixing. This operator (see Section 4) opens an invisible decay channel for the Higgs boson into an HNL and a light neutrino, with a rate given by where we have neglected the light neutrino mass. Note that no flavour index has been specified; this decay affects all three flavour copies of the operator, so the squared Wilson coefficient in the equation above stands for the sum of all three squared coefficients: Neutral bosonic current. A vertex between a Z boson and two HNLs is introduced (Section 5.1), thus mediating invisible decays for the Z and for neutral mesons, both pseudoscalar, P 0 , and vector, V 0 . The corresponding rates are Here, M Z stands for the Z mass, G F for Fermi's constant, and x N is the ratio of the HNL mass with respect to its parent particle mass. M P and f P are the neutral pseudoscalar meson mass and decay constant, 7 and similarly for M V and f V for the neutral vector mesons. g V accounts for the coupling of the Z boson to the quark content of the meson under consideration. In particular, for the Υ(1S), g Υ = √ 2 1 2 − 2 3 s 2 W , where s W is the sine of the weak mixing angle (see Ref. [140] for more details).
Charged bosonic current. This operator mediates the most relevant processes for HNL production, mainly in charged meson decays. Furthermore, there is a direct relation between the Wilson coefficient and the effective mixing they induce (Eq. (5.3)). Thus, the usual expressions for decay rates involving HNLs apply for this effective operator. We refer the reader to Ref. [140] for a list of the relevant expressions.
However, the absence of neutral currents affects some HNL decay channels, which, in the case of standard mixing, involve NC-and CC-like diagrams. This is the case for the N → νee and N → νµµ processes. The corresponding rate mediated by this operator is analogous to that of the decay N → νqq′ in the case of standard mixing [141], and reads with x α ≡ m ℓα /M N and λ(x, y, z) is the Källen function.
Neutral four-fermion interactions. The operators O uu , O dd (Eq. (7.1)) and O QN (Eq. (7.3)) mediate invisible decays of neutral mesons into two HNLs. Although the chiralities of the quark fields are different, the expressions for the widths will be the same, since the vector contribution vanishes. They are given by The couplings g P and g V account for the overlap of the quark content of the involved meson and the quark structure of the operator under consideration. They are summarized in Table 2 for the most relevant mesons. We refer the reader to Ref. [140] for the values of most of the decay constants, as well as for the angles θ 0 and θ 8 , that parametrize the mixing in the η sector. For the heavier mesons we employ f J/ψ = 1.8 GeV 2 [142] and f Υ = 9.2 GeV 2 [143]. A four-lepton interaction is induced by the operator in Eq. (7.4), inducing µ and τ decay into lighter leptons. The corresponding width is given by: Charged four-fermion interactions. The operators in Eqs. (7.5), (7.7) and (7.11) induce charged meson decays into HNLs and vice versa (we only display here the pseudoscalar case, as it is the main HNL production source). The operator O duNℓ provides the same rates as in the standard mixing case (see Ref. [140]) but controlled by the "effective mixing" where i, j are the flavours of the quarks composing the meson, and V ij is the corresponding CKM element. The operators O LNQd and O QuNL exhibit the same decay widths for mesons into HNLs and vice versa: where m i and m j are the masses of the quarks that compose the meson under consideration.