GeV-scale neutrinos: interactions with mesons and DUNE sensitivity

The simplest extension of the SM to account for the observed neutrino masses and mixings is the addition of at least two singlet fermions (or right-handed neutrinos). If their masses lie at or below the GeV scale, such new fermions would be produced in meson decays. Similarly, provided they are sufficiently heavy, their decay channels may involve mesons in the final state. Although the couplings between mesons and heavy neutrinos have been computed previously, significant discrepancies can be found in the literature. The aim of this paper is to clarify such discrepancies and provide consistent expressions for all relevant effective operators involving mesons with masses up to 2 GeV. Moreover, the effective Lagrangians obtained for both the Dirac and Majorana scenarios are made publicly available as FeynRules models so that fully differential event distributions can be easily simulated. As an application of our setup, we numerically compute the expected sensitivity of the DUNE near detector to these heavy neutral leptons.


Introduction
The evidence for neutrino masses and mixings from the neutrino oscillation phenomenon demands an extension of the Standard Model (SM) of particle physics so as to accommodate the experimental results. Arguably, the simplest of such extensions is to add fermion singlets to the SM particle content. Indeed, the inclusion of these right-handed neutrinos would make the neutrino sector equivalent to its charged-lepton counterpart and allow for neutrino Yukawa couplings in complete analogy to the other fermions of the SM. However, being complete singlets of the SM gauge group, the novel and distinct option of a Majorana mass term is also open for them. This Majorana mass term would not only include a new source of particle number violation, possibly related to the origin of the observed baryon asymmetry of the Universe (BAU), but also include a new energy scale in the Lagrangian, not connected with electroweak symmetry breaking. As such, there is no solid theoretical guideline for the value of this new physics scale.
An attractive possibility is that the Majorana mass scale is much larger than the electroweak scale, possibly close to the Grand Unification scale, leading to the celebrated type-I Seesaw mechanism [1][2][3][4]. Its most appealing feature is that the smallness of neutrino masses is very naturally explained even with order one Yukawa couplings, since it is inversely proportional to the large Majorana mass of the right-handed neutrinos. Furthermore, the Seesaw mechanism, unlike the SM, is also able to account for the observed BAU via leptogenesis [5]. Nevertheless, while a very high Majorana mass scale can very naturally accommodate the extreme lightness of neutrino masses, its presence would significantly destabilize the Higgs mass, worsening the Higgs hierarchy problem [6,7]. Thus, while naturalness arguments favour large Majorana masses to explain the light neutrino masses, lighter scales are instead preferred to accommodate the observed Higgs mass.
Other variants of the original Seesaw mechanism, such as the inverse [8][9][10] or linear [11] Seesaws, naturally explain the lightness of neutrino masses through an approximate lepton number symmetry [12][13][14] instead. Thus, they can be realized at lower energy scales without introducing a Higgs hierarchy problem. Low-scale versions of the leptogenesis mechanism are also found to successfully account for the observed BAU [15][16][17]. Thus, the phenomenology associated to all possible options for the Majorana mass scale should be investigated and compared with experimental observations so as to probe the new physics underlying the observation of neutrino masses and mixings.
The main consequence of lowering the Majorana mass scale in a Seesaw mechanism is that mainly-sterile neutrinos or "heavy neutral leptons" (HNLs) appear in the particle spectrum. If sufficiently light, these HNLs will be kinematically accessible to experiments and can thus be produced and searched for. Given the singlet nature of the right-handed neutrinos, their only interactions are the weak ones, inherited from their left-handed neutrino counterparts via mixing. Thus, the mixing of the HNLs with the electron, muon and tau neutrinos can be probed and constrained as a function of the HNL mass by searching for their production and decay in association with the corresponding charged leptons. These searches range from studying their impact in neutrino oscillations [18][19][20][21][22][23][24][25][26] when they are too light to decay visibly, to collider signals [27][28][29][30][31][32][33][34][35][36][37][38] for the highest accessible HNL masses. For even higher masses, their mixing can still be constrained through deviations of unitarity of the PMNS matrix in flavor and electroweak precision observables [39][40][41][42][43][44][45][46][47][48].
For intermediate HNL masses M , between the MeV and GeV scales, searches at beam dump experiments or near detectors of neutrino oscillation facilities, where they can be produced via meson decays and detected through their visible decays, can set very stringent constraints [14,28,[49][50][51][52][53][54][55][56][57][58][59][60][61][62]. Indeed, current bounds are even getting near the expectation for the "vanilla" type-I Seesaw without a lepton number symmetry protection of the light neutrino masses m ν , where the mixing scales as θ 2 ∼ m ν /M . Furthermore, the masses and mixings leading to successful generation of the BAU via low-scale leptogenesis are also accessible through these searches [63][64][65]. In this regime, both the production and decay of the HNL depend crucially on its interactions with mesons. While these have been studied previously, significant discrepancies can be found in the literature [28,49,53,66] in the branching ratios of the relevant channels. The aim of this work is to clarify such discrepancies and provide a tool for these important searches. With that goal in mind, we derive the effective theory description of the HNL interactions, with particular emphasis on the effective operators involving mesons, which control HNL production and decay via leptonic and semileptonic processes. We do this both for a Majorana HNL as well as for the Dirac scenario, motivated by the inverse and linear Seesaw variants. Furthermore, our results have been collected in two FeynRules [67] models that have been made publicly available (see ancillary files) so that not only the total branching ratios can be computed, but also differential event distributions can be easily simulated by interfacing the output of FeynRules with event generators such as MadGraph5 [68]. Finally, while the present work focuses on the low-energy theory, our FeynRules implementation is more general and includes an option to replace all mesons with quarks, so they may also be used to study HNL phenomenology in collider searches at higher energies.
As an application of our framework, we compute the expected flux of HNLs at the proposed DUNE [53,54] near detector, and compare our full numerical simulation with the approximation of rescaling the massless neutrino fluxes. A significant enhancement due to the larger boost in the beam direction for the HNLs is found. We also compute the expected number of decays inside the DUNE near detector into several decay channels, and use that to estimate the sensitivity of DUNE to the HNL mixing with the charged leptons. This paper is organized as follows. In Section 2 we introduce the Seesaw Lagrangians, both in the Majorana and Dirac cases, and review the weak interactions that the HNL will inherit from the left-handed neutrinos via mixing. In Section 3 we concentrate on the meson interactions and derive all the relevant effective operators containing HNLs. In Sections 4 and 5 we summarize all the relevant production channels and subsequent decays of the HNLs. In Section 6 we present our results for the expected HNL fluxes at the DUNE near detector together with an estimate of its sensitivity. Finally, in Section 7 we draw our conclusions and summarize the results.

The full Lagrangian of the theory at high energies
Once the SM is extended with n extra right-handed neutrinos N R , Lorentz and gauge invariance allow the inclusion of Yukawa couplings to the lepton doublets (Y ν ) as well as Majorana masses for the heavy singlets (M ). In the basis where the Majorana mass terms are diagonal, the corresponding Lagrangian reads: where L L,α stands for the SM left-handed lepton doublet of flavor α, ϕ is the Higgs field, ϕ = iσ 2 ϕ * and N c R,j ≡ CN t R,j , with C = iγ 0 γ 2 in the Weyl representation we adopt. Once the Higgs develops its vacuum expectation value v/ √ 2 upon electroweak (EW) symmetry breaking, the full neutrino mass matrix in the basis (ν L , N c R ) can be written in blocks as: The full unitary rotation U that diagonalizes the mass matrix will have dimensions (3 + n) × (3 + n). Neutrino masses are obtained upon diagonalization, as well as the mixing between the active SM neutrinos and the new heavy states introduced. In particular, the spectrum is composed of 3 light "SM-like" neutrino mass eigenstates (ν i ), and n heavier and mostly sterile neutrinos (N i ). Alternatively, and motivated by low-energy Seesaw realizations such as the inverse [8][9][10] or linear [11] versions, we will also consider the case in which the extra sterile neutrinos have Dirac (or pseudo-Dirac) masses. In these scenarios, 2n extra singlets are added in Dirac pairs N L,j , N R,j (j = 1, . . . n). Neglecting the small lepton-number violating terms (that would eventually source the light neutrino masses), we are left with the following Lagrangian: In this case, the mass matrix in the basis (ν L , N c R , N L ) would be given by: Regardless of the Dirac or Majorana character of the heavy neutrinos, the flavor states will thus correspond to a combination of the light and heavy states: where we have introduced the mass eigenbasis n = (ν, N ) with index i that runs over the light and heavy mass eigenstates. The leptonic part of the electroweak Lagrangian can be written as: where P L and P R are respectively the left and right projectors, c w ≡ cos θ w , s w ≡ sin θ w (θ w being the SM weak mixing angle), and The heavy neutrinos can also interact with the quark sector through the charged and neutral current interactions. Thus, the corresponding weak interactions between quarks are reviewed below for convenience: and Here, Q q and T q 3 stand for the electric charge and the isospin of quark q in the interaction vertex (from now on the index q will be dropped for simplicity), and V qq ′ is the corresponding element of the CKM mixing matrix.
Finally, for the derivation of the effective theory in Sec. 3 it will be useful to separate both currents in their vector and axial parts. This way, the Z current can be decomposed as Analogously, the W current may be written as where its vector and axial parts are given by

Effective low-energy Lagrangian including mesons
In order to compute the production of the heavy neutrinos through meson decays, as well as neutrino decays to lighter mesons, we need to introduce effective interactions between the neutrino and meson fields. In this section we derive such interactions, integrating out the W and Z bosons and introducing the relevant meson decay constants and hadronic matrix elements. We compute the amplitudes for low-energy processes involving these vertices, so as to extract the corresponding effective operators. Moreover, FeynRules [67] models with these effective interactions have been made publicly available (as ancillary files to this work), making possible the generation of fully differential event distributions. Note that, while the formalism used in this section is applicable to any number of extra heavy states, only one heavy neutrino has been included in the FeynRules model files for the sake of simplicity. Although the introduction of just one heavy neutrino cannot explain the measured neutrino masses and mixing parameters, such simplified models are useful to study the phenomenology of HNLs, since it will be dominated by the lightest of the extra states.
As a first step, we review the relevant decay constants and matrix elements, and introduce our notation. Throughout this section, the formalism we use is suitable for mesons with masses up to approximately 1 GeV. However, leptonic and semileptonic decays of heavier charmed mesons can constitute a dominant contribution for heavy neutrino production, depending on its mass. Such processes will also be considered here, for the channels with a significant branching ratio into neutrinos (e.g., D s → N ℓ). For even higher neutrino masses (produced typically at collider experiments), a perturbative description of the neutrino decay into quark-antiquark pairs (with subsequent hadronization) would be more suitable.
We adopt a definition of the meson decay constants such that f π = 130 MeV, namely: for the pseudoscalar (P ) and vector (V ) mesons respectively. Here, p µ stands for the momentum of the pseudoscalar meson and ϵ µ for the polarization of the vector meson (note that, with this definition, the decay constants f V have units of [E] 2 ). The corresponding currents j A a,µ , j V a,µ are defined as: (3.5) In this notation, the set {λ a } corresponds to linear combinations of the eight Gell-Mann matrices (generators of SU (3)) plus the identity, normalized such that For convenience, explicit expressions for the generators are provided in Appendix A, while the decay constants most relevant for the effective couplings considered in this work are summarized in Tab. 1.

Pseudoscalars Vectors
Decay constants Rotation angles  [69], while those for vector mesons have been computed as described in Appendix C. Right. Decay constants for the η 0 and η 8 , and angles that parametrize the rotation to the physical basis, taken from Ref. [70] (see text for details). Note that in Ref. [70] the authors use a different normalization for the current definitions than the one adopted in this work. However, this does not affect our result since they provide their results in terms of the ratios f 8 /f π and f 0 /f π , which remain unaffected by an overall normalization factor.

Pseudoscalar mesons
3.1.1 Neutral mesons: π 0 , η, η ′ The quark content of the neutral pseudoscalar mesons will correspond to linear combinations of the diagonal generators λ 0 , λ 3 and λ 8 . Substituting the explicit expressions for the generators into Eq. (3.3) we obtain: The neutral pion can be directly identified with the current j A 3,µ , being the neutral member of the SU(2) triplet of pseudo-Goldstone bosons from the flavor symmetry between upand down-quarks. Conversely, the η and η ′ mainly correspond to the currents j A 8,µ and j A 0,µ respectively, although with significant mixing among them, as discussed in detail below. These neutral mesons can be produced or decay through neutral current interactions mediated by the Z boson. Thus, in order to obtain their effective interactions with neutrinos we start from the Fermi theory after integrating out the Z, inserting the decay constant of the corresponding meson. The Z axial current in Eq. (2.13) can be expressed as a linear combination of the neutral axial currents as: At low energies, the amplitude, for example, for π 0 → n inj would read: whereū i and v j are the corresponding spinors for the neutrino mass eigenstates. Substituting the Z current from Eq. (3.8) and the corresponding hadronic matrix element from Eq. (3.1), and introducing Fermi's constant, the amplitude is given by: where p µ is the 4-momentum carried by the pion. Translating the momentum into a derivative, we can write down, in configuration space, the effective operator that leads to the amplitude in Eq. (3.11): Furthermore, if all particles are on-shell, it is possible to apply Dirac's equation to obtain Yukawa couplings proportional to the neutrino masses: Since the coupling is proportional to the masses of the neutrinos, the coupling to the heavy states will dominate the interaction, in complete analogy to the chiral enhancement of the charged pion decay π → µν µ versus π → eν e . Similarly, the operators associated to the other neutral pseudoscalar currents, for onshell particles, can be obtained as: However, unlike in the π 0 case, the η and η ′ mesons mix significantly and do not correspond exactly with the quark content of the η 8 and η 0 defined through the corresponding currents in Eq. (3.7). Thus, a change of basis must be performed in order to obtain the effective vertices for the physical states. We adopt the usual parametrization for this change of basis, with two angles, θ 0 and θ 8 (see e.g. Ref. [70]), and define: The values for f 0 , f 8 , θ 0 and θ 8 have been taken from Ref. [70] and are summarized in Tab. 1 for convenience. Through this change of basis, the currents for the η and η ′ can be obtained as combinations of the j A 0,µ , j A 8,µ currents as Therefore, the relevant operators in the mass basis will read The normalized combinations of generators that reproduce the quark content of the π ± and K ± are: Thus, from Eq. (2.16) we get that The amplitude for π − → ℓ −n is obtained after integrating out the W boson, following the same procedure used to derive the effective vertex for the π 0 →nn decay in the previous section: After introducing the W current defined in Eq. (3.23) and evaluating the hadronic matrix element, the amplitude reads: In the same fashion as before, we translate this amplitude to an effective operator in configuration space, with the 4-momentum p µ as a derivative acting on the leptonic current: Once again, if all the particles involved are on-shell, it is possible to obtain Yukawa couplings proportional to the fermion masses via Dirac's equation: This procedure can be repeated for the charged kaons, obtaining the same result once the corresponding decay constant and CKM element are introduced: So far we have restricted ourselves to mesons which contain only the three lightest quark flavors. Nevertheless, these results can be generalized to the D ± and D ± s mesons. The corresponding effective operators read:

Vector mesons
3.2.1 Neutral mesons: ρ, ω, ϕ As for the pseudoscalar case, the vector currents associated to the generators can be expressed in terms of the u, d and s quarks as Considering their respective quark contents, the corresponding normalized currents for the ρ 0 , ω and ϕ mesons are given by: The production and decay of the vector mesons take place via the vector component of the Z current, Eq. (2.12), which can be written as the following linear combination of the vector meson currents: After integrating out the Z boson, the amplitude for the ρ 0 →nn process reads: Introducing the vector Z current defined in Eq. (3.33) and evaluating the matrix element according to Eq. (3.2), we get: where ϵ µ is the polarization vector of the ρ 0 meson. It is then immediate to extract the effective operator in configuration space: Analogously, for the other two neutral vector mesons we obtain: In complete analogy to the charged pseudoscalars, the charged vector meson currents are given by: 40) and the vector component of the W current from Eq. (2.15) can be written as: The computation of the effective operators is done exactly in the same way as for the charged pseudoscalar case. The amplitude for the ρ − →nℓ − process reads: Thus, we finally obtain and, equivalently, for the K * ,± meson we get

Semileptonic meson decays
Some mesons exhibit non-negligible branching ratios for semileptonic decay channels into neutrinos, charged leptons and lighter mesons. These can even dominate over the two-body leptonic decays if the mass of the heavy neutrino is not large enough to sufficiently enhance the latter, and thus must be taken into account. After integrating out the W boson, the amplitude for the P → Dnℓ decay (where P and D stand for generic parent and daughter mesons, respectively) reads: where j V W,µ is defined in Eq. (2.15). This hadronic matrix element is usually expressed in terms of two form factors, f + and f − [71]: where V qq ′ is the CKM element corresponding to the quarks which interact with the W in the hadronic transition, while p µ ≡ p D µ + p P µ is the sum of the 4-momenta of the parent and daughter mesons and q µ ≡ p D µ − p P µ is the 4-momentum transfer between them. Thus, the amplitude can be written as: In what follows, it becomes convenient to express this in terms of the 4-momenta of the daughter meson, p D µ , and of the leptonic pair, p nℓ µ : Note that we have not specified the electric charges of the involved mesons. In fact, this amplitude describes all the processes allowed by charge conservation (P − → D 0n ℓ − and P 0 → D +n ℓ − , as well as their CP-conjugates). However, it should be stressed that, even though electromagnetic contributions to these hadronic form factors are generally small, in some cases the numerical parameters they contain might be slightly different depending on the charge of the mesons, since they come from fits to different datasets. From this amplitude it is possible to extract the corresponding effective operator in configuration space, writing the 4-momenta as derivatives: where ϕ P and ϕ D are the parent and daughter meson fields, respectively. Once more, if the involved fields are on-shell, it is possible to apply Dirac's equation and substitute the derivative acting on the leptonic current by terms proportional to their masses. The resulting operator reads:

Form factors
Many parametrizations for the hadronic form factors are available in the literature, most of which are given in terms of f + and f 0 . The former was defined together with f − in Eq. (3.46), while the latter can be related to f + and f − via The semileptonic decays we will be mostly interested in are K → πnℓ and D → Knℓ. For the former we employ a linear parametrization, as in Ref. [72], according to which Conversely, in the case of the D → Knℓ decay we make use of a "pole" parametrization [71]: The values used for the form factor parameters are summarized in Tabs. 2 and 3.
0.7647 −0.066 −2.084 Table 2. Parameters entering our form factor definitions for the semileptonic D → Knℓ decays. Table 3. Parameters entering our form factor definitions for the semileptonic K → πnℓ decays. Note that we make use of different parameters for the decays of charged and neutral kaons, following Ref. [69].
We have included these form factors in our FeynRules model, and numerically checked with MadGraph5 [68] that our implementation reaches an agreement of at least a 95% with the measured branching ratios for the SM decay channels K → πνℓ and D → Kνℓ [69]. For convenience we provide two separate implementations for such couplings, as explained in detail in Appendix D.

Production of Heavy Neutral Leptons from meson decays
In this section we provide the expressions for the production of a heavy neutrino N 4 of mass M 4 via meson decays. We have computed them employing the Feynman rules derived from the effective operators obtained in Sec. 3, and verified their agreement with the simulations generated via MadGraph5 using the implementation of our model in FeynRules. In order to do so, we have diagonalized explicitly the full mass matrix and expressed the ensuing mixing matrix in terms of the original Yukawa couplings. Further details on this diagonalization can be found in Appendix B.

Two-body leptonic decays
The generic expression for the leptonic decay of a charged pseudoscalar meson P of mass m P is given by [28,49,66,73,74] where the values of f P are given in Tab. 1, and we have defined

Three-body semileptonic decays
The decay width for the semileptonic decay of a parent pseudoscalar meson P into a daughter pseudoscalar D, a charged lepton ℓ α and a heavy neutrino N 4 is given by [49,66,73,74] where C D = 1 in all cases under consideration, except for K ± → π 0 N 4 ℓ ± α , for which C D = 1 √ 2 . The integrals I P D i are expressed in terms of the form factors f P D + (q 2 ) and f P D 0 (q 2 ) defined in Sec. 3.3.1.
5 Decays of Heavy Neutral Leptons into SM particles

Two-body decays
Here we provide general expressions for the decay widths of a heavy neutrino N 4 of mass M 4 into final states including pseudoscalar and vector mesons separately. We have computed them employing the Feynman rules derived from the effective operators obtained in Sec. 3, and verified their agreement with the simulations generated via MadGraph5 using our model implementation in FeynRules. Throughout this section, we will neglect the masses of the light neutrinos for simplicity.

Pseudoscalar mesons
The generic expression for the heavy neutrino decay width into a neutral pseudoscalar meson P is given by where we have defined x P ≡ m P /M 4 , and according to the parametrization used to describe the η − η ′ mixing in Sec. 3.1.1. Using the parameters provided in Tab. 1, this leads to the "effective decay constants" f η ≃ 81.6 MeV and f η ′ ≃ −94.6 MeV. Finally, note that the sum j in Eq. (5.1) runs over the three light neutrino mass eigenstates, since they cannot be individually identified. However, at leading order in U α4 , this is equivalent to a sum running over the three active flavors, since On the other hand, the decay width into a charged pseudoscalar meson P ± is given by

Vector mesons
In the case of neutral vector mesons, the decay width reads: where the decay constants f V are again summarized in Tab. 1.

Three-body decays
Heavy neutrinos may also decay into three body final states either purely leptonically or semileptonically. The latter include N 4 → π + π 0 ℓ − , N 4 → π 0 π 0 ν and N 4 → K + π 0 ℓ − . However, their respective contributions are dominated by N 4 → ρ + ℓ − , N 4 → ρ 0 ν and N 4 → K * ,+ ℓ − respectively, already included in the previous section. This can be seen from the data from τ decays, since the hadronic matrix elements involved in the semileptonic decays would be the same. Indeed, the branching ratio of τ − → νπ − π 0 is 25.49%, while the contribution which does not correspond to τ − → νρ − is negligible: (3.0 ± 3.2) · 10 −3 [69]. We will thus review here only the three-body purely leptonic decays N 4 → ℓℓν and N 4 → ννν, taken from Refs. [28,49,66,73]. The invisible decay of the heavy neutrino reads [28,49,66,73] where we have summed over all possible light neutrinos in the final state. For the three-body decays involving charged leptons in the final state, we will distinguish between two cases. If the heavy neutrino decays into two leptons of the same flavor β, there are both W and Z mediated diagrams contributing to the amplitude. The total decay width can be expressed as [28,49,66,73] and we have defined the functions On the other hand, the decay of the heavy neutrino into two leptons of different flavor is only mediated by the W interaction. In the limit in which one of the charged lepton masses can be neglected, the corresponding decay width simplifies to [49,66] Note that this expression corresponds to a Dirac neutrino decay; for Majorana neutrinos there would be a second contribution proportional to |U β4 | 2 since there are two diagrams allowed, each of them proportional to a different mixing matrix element.

Decays to 4 or more bodies
Finally, for HNL masses above 1 GeV, the appropriate description of the hadronic final states transitions from the effective theory with the different meson resonances included in the previous sections, to quark production in the final state with subsequent hadronization, more suitable for perturbative QCD. For reference, the τ − , with its 1.78 GeV mass, is precisely at the transition region. Indeed, it shows a 10.8% and 25.5% branching ratio to the ν τ π − and ν τ ρ − channels respectively. But also a 9.3% branching ratio to both ν τ π − 2π 0 and to ν τ 2π − π + and even 4.6% and 1.0% branching ratios to ν τ 2π − π + π 0 and ν τ π − 3π 0 respectively [69]. These last decay modes, with three or more mesons in the final state, are more suitably described from the underlying quark interactions with a subsequent correction to account for the hadronization process: 1 + ∆ QCD ≡ Γ (τ → ν τ + hadrons) Γ tree τ → ν τ + u +d + Γ tree (τ → ν τ + u +s) (5.14) with [75] ∆ QCD = α s π We adopt the same approach as Ref. [66] (see also [59,76]) and use Eq. (5.15) to account for the hadronization of the HNL decays N 4 → ℓ α ud and N 4 → ℓ α us for HNL masses above 1 GeV. We also apply the same correction to the neutral current decays N 4 → νqq with q = u, d, s. However, we add a phase space suppression factor 1 − 4m 2 K /M 2 4 for the N 4 → νss channel since it would otherwise overestimate its importance for M 4 = 1 GeV, where the phase space prevents two K in the final state. For the running of α s we follow the dedicated review in Ref. [69]. The difference between these fully inclusive hadronic final states and the HNL decays to specific mesons discussed above will provide an estimate of the HNL decays to 3 or more mesons. We have tested this procedure for the τ decays and reached good agreement with its tabulated branching ratios. Figures 1 and 2 show the branching ratios for the different decay channels of the heavy neutrino, as a function of its mass, for two different cases: degenerate mixings to all lepton flavors (|U e4 | 2 = |U µ4 | 2 = |U τ 4 | 2 ), and in the case when only one of the mixing matrix elements is non-zero. The labels ℓ ∓ hadr. and νhadr. stand for N 4 decays, mediated by charged and neutral currents respectively, with 3 or more mesons in the final state.

Discrepancies with previous literature
The decay widths of a HNL into mesons, neutrinos and leptons have been derived several times in previous literature; for an incomplete list see e.g. Refs. [28,39,49,53,66,73]. Here we summarize the main discrepancies and differences found between our results and some of these works: 1. Overall, we find a relatively good agreement with Ref. [66] for the meson decay constants and for most vertices involving heavy neutrinos, with the exception of the couplings to ω and ϕ mesons, for which we find different expressions in terms of sin 2 θ w (see our Tab. 4, in comparison with Tab. 9 in Ref. [66]).
2. We find that the expressions in Ref. [49] for the HNL decay into vector mesons have an extra factor 2 with respect to our results, both for the neutral and the charged channels. Also, their expressions for the decay into neutral vector mesons seem not to include a dependence on sin 2 θ w (see our definitions for g V in Tab. 4). Finally, there are significant differences in the values reported in Ref. [49] for the neutral pseudoscalar meson decay constants f η and f η ′ .
3. We find that the expressions for the HNL decay into a light neutrino and a neutral pseudoscalar meson in Refs. [28,73] have an extra factor 2 in the denominator, as the authors of Ref. [66] pointed out. Branching ratios for the heavy N 4 as a function of its mass, obtained under the assumption that only its mixing to one lepton flavor is non-zero, as indicated by the labels in each row. Left (right) panels correspond to decays without (with) light neutrinos in the final state. The decay channels into semileptonic final states are not shown as their branching ratio is expected to be negligible in the range of masses considered here. 4. We also find significant discrepancies with the HNL decays to neutral vector mesons in Ref. [28], which were also already pointed out in Ref. [66].
5. Regarding Ref. [53], which was published more recently, we again find some discrepancies on the branching ratios for vector mesons: our results show a significantly higher branching ratio for the decay channels N 4 → ρℓ, and lower branching ratios for the N 4 → ων and N 4 → ϕν decays (as can be seen from the comparison between our Fig. 1 and their Fig. 1). These discrepancies can be partially explained by the different couplings we obtain for the neutral vector meson couplings (see our Tab. 4, in comparison with the values given below Eq. (3.12) in Ref. [53]) and possibly by the different values used for the corresponding decay constants. Indeed, there are also significant discrepancies in the literature for the choices of the vector meson decay constants. We thus clarify our choice in Appendix C.

Heavy Neutral Leptons at DUNE
In the remainder of this work, we use the effective theory derived in the previous sections to compute the expected heavy neutrino flux at the DUNE near detector (ND), as well as the expected number of decays for different channels. From now on we will assume a Dirac HNL; in the Majorana case, the heavy neutrino decay widths, and thus the number of events, would increase in a factor of 2. In all our calculations, we consider a ND geometry as described in the DUNE Technical Design Report (TDR) [77]. The ND complex will be located 574 m downstream from the neutrino beam source, and will include three primary detector components: a liquid Argon Time Projection Chamber (LArTPC) called ArgonCube; a high-pressure gaseous TPC surrounded by an electromagnetic calorimeter (ECAL) in a 0.5 T magnetic field, called the Multi-Purpose Detector (MPD); and an onaxis beam monitor called System for on-Axis Neutrino Detection (SAND). We will consider the detector volume corresponding to the MPD 1 for which the beam-induced background is smaller given its lower density. All calculations presented in this section have been performed using the nominal beam configuration and luminosity envisioned for DUNE [77]: 120 GeV protons and 1.1 · 10 21 protons on target (PoT) per year, divided equally into positive and negative horn focusing modes, which yields a total of 7.7 · 10 21 PoT over 7 years of data taking, which is expected to start in 2027. The simulation of the meson production in the target has been done as follows. For pions and kaons, we use the results of the detailed GEANT4 [80][81][82] based simulation (G4LBNF) of the LBNF beamline developed by the DUNE collaboration [77]. The simulation includes a detailed description of the geometry, including the 1.5 m long target, three focusing horns, decay region, and surrounding shielding. The DUNE collaboration provides both neutrino and antineutrino mode predictions, generated for a 120 GeV primary proton beam. For positive horn focusing mode (PHF) we use the results of the full simulation to calculate the predicted event rate at the DUNE ND, while for negative horn focusing (NHF) mode we scale the event rates from PHF mode based on the flux ratios between π − /π + and K − /K + as predicted by G4LBNF.
However, G4LBNF does not include the production of D, D s and τ leptons. Thus, in this case Pythia (v 8.2.44) [83] was used to create a pool of events and predict production rates for proton collisions at various momenta, and a GEANT4-based simulation was subsequently used to predict proton inelastic interactions with 120 GeV primary protons impinging on the target. For each inelastic interaction, we randomly pick a Pythia event from the pool of events generated at the corresponding momentum, with a weight proportional to the rate predicted by Pythia. In doing this, we neglect the effect of the magnetic horns since these heavy particles decay very promptly and, therefore, it is safe to assume that their production will be similar for the PHF and NHF modes. The average number of parent mesons and τ leptons per PoT produced in the target 2 are listed in Tab. 5.  Table 5.
Average number of positive and negative parent mesons and τ leptons P per PoT produced in the target. Tab. 6 summarizes the different HNL production channels that have been included in our analysis. This set contains the dominant leptonic and semileptonic decays into heavy neutrinos of the parent mesons (π, K, D and D s ) produced in the target. Moreover, since the D and D s decay very promptly and have sizable branching ratios to τ leptons, a significant τ production rate is expected. This provides an additional production mechanism for HNL masses below the τ mass controlled by |U τ 4 | 2 , allowing DUNE to significantly improve the sensitivity to this more elusive mixing matrix element. All decay modes of the τ could allow to produce a HNL in the final state, provided that it is kinematically allowed. Nevertheless, we have opted to conservatively consider only the τ decay modes Unlike for the production from meson decays, we did not provide their explicit expressions in Sec. 4, since they would be the same as for the corresponding N 4 decays in Eqs. (5.6), (5.4) and (5.13), respectively. The decays with 3 or more mesons in the final state have been neglected since the phase space is reduced for the production of a massive particle and the simulation of the HNL kinematics is more challenging for these channels (see Sec. 5.3).
Once we have obtained an expected flux of HNL entering the detector, this is then matched to the 22 different decay modes into SM particles studied in Sec. 5 (and shown in Figs. 1 and 2) according to their corresponding branching ratios to obtain the expected signal at the detector.
In the remainder of this section we first illustrate the impact on the detector acceptance due to the boost of the HNL, and then we compute the expected number of heavy neutrino decays inside the DUNE ND to estimate its sensitivity.

The effect of the HNL mass on the detector acceptance
The effect of the boost in the beam direction is more efficient for particles with smaller velocities. Therefore, a larger detector acceptance is obtained when the effect of the heavy neutrino mass is properly included in the flux simulations, compared to estimations of the HNL flux based on the massless neutrino distributions. In order to illustrate the effect of Parent 2-body decay 3-body decay Table 6. List of 2-body and 3-body decays into HNLs, for the parent particles considered in this work. The decay channel τ − → π − π 0 N 4 is simulated via the approximation τ − → ρ − N 4 , ρ − → π − π 0 as discussed in the text.  this figure, the increase in acceptance is considerable: up to a factor of two for 200 MeV neutrinos from kaon decays, and up to a factor of three for 1 GeV neutrinos coming from D decays. The effect of the boost will also lead to a distortion in the expected spectra due to the different dependence of the detector acceptance with the neutrino energy, which can be seen from the comparison of the shape of the light and dark histograms in each panel. The net result is a relative increase in the number of neutrinos at low energies that enter the ND, given their smaller velocities and hence stronger collimation.

Parent 2-body decay 3-body decay
Finally, Fig. 4 shows the total detector acceptance, after integrating over the neutrino energy, as a function of the HNL mass. Note that the acceptance is expected to be different depending on the parent meson that produced the neutrino due to the effect of the horns: while pions are typically very well-focused at a long-baseline experiment, this is not the case for heavier mesons, which are not only harder to focus due to their larger masses but also decay much faster. This effect is most significant for D and D s mesons, which decay very promptly and therefore are practically unaffected by the horn focusing system. In Fig. 4 we show different lines for neutrinos obtained from different meson decays, as indicated by the labels. As can be seen, as the heavy neutrino mass approaches the production threshold and its velocity decreases accordingly, the acceptance grows very rapidly given the stronger boost in the beam direction. Notice however that also the phase space is decreasing and hence the number of total HNL events will also be reduced.

Expected sensitivity to HNL decays
Once the flux of heavy neutrinos dϕ N /dE N that reach the ND has been computed numerically as a function of the neutrino energy E N , the total number of expected neutrino decays into a given decay channel c inside the DUNE ND can be expressed as where BR c is the branching ratio of the corresponding decay channel and P (E N ) stands for the probability of the heavy neutrino decaying inside the ND (which depends on the boost factor and therefore on the neutrino energy): Here, Γ is the total decay width of the heavy neutrino in its rest frame, while γ = E N /M 4 , β = | ⃗ p N |/E N and ⃗ p N stands for the neutrino momentum. L is the distance between the HNL production and the ND while ∆ℓ det is the length of the HNL trajectory inside the detector. From Eq. (6.2) it is easy to see that the neutrino must be sufficiently long-lived to reach the ND, or otherwise the number of decays will be exponentially suppressed. This will be the case for large enough energies and small enough matrix elements U α4 , which correspond to the most interesting region of the parameter space. In this limit, ΓL ≪ γβ, and the decay probability can be further approximated as Nevertheless, this approximation does not hold anymore for large masses and mixings, and thus we will employ Eq. (6.2) to compute the probability of the HNL decaying inside the detector.
Given that the neutrino flux entering the detector will be directly proportional to its aperture, and that the probability for the neutrino to decay inside the ND is proportional to ∆ℓ det , it is easy to see that the sensitivity to heavy neutrino decays will scale with the volume of the ND. In order to compute the final sensitivity to HNLs, a detector simulation should be performed, including relevant background contributions from SM neutrino interactions in the ND. Such a fully detailed detector simulation is beyond the scope of this work, where we rather show an estimate to the sensitivity of DUNE as an application of the methods derived in the previous sections. The main source of background for this search comes from neutrino interactions in the detector volume, and is very significant. In Ref. [53] the background rates for Argon were estimated at ∼ 3 · 10 5 events/ton/10 20 PoT. Fortunately, SM neutrino events present a very different topology than that of heavy neutrino decays, and a series of kinematic cuts can heavily reduce the expected background and bring it down to a negligible level. This was the case, for example, for the T2K near detector HNL search performed in Ref. [84] (which also used a gas TPC). Therefore, following Refs. [53,84], hereafter we will assume that this is achievable and show our expected sensitivity contours to heavy neutrino decays under the assumption of no background. We also assume that the cuts applied to reduce the background will translate into similar signal efficiencies in our case as those obtained in Ref. [84]. Although the efficiency will eventually depend on the mass of the HNL and the considered decay channel (see Fig. 4 in Ref. [84]), here we use 20% as an educated guess. Finally, for our sensitivity contours we estimate the 90 % confidence level (CL) sensitivity on the signal following the Feldman and Cousins [85] prescription for a Poisson distribution with no background and under the hypothesis of no events being observed, which corresponds to the expected number of signal events being smaller than 2.44.
Before showing our sensitivity contours, we show in Fig. 5 an example to illustrate the relative importance of the different HNL production mechanisms on the results. In this example, we show the contours obtained under the assumption that the heavy neutrino mixes primarily with the e sector. The different regions show the contributions obtained when the heavy neutrinos are produced from the decays of a given parent meson, as indicated by the labels. The signature in this case would be electron-positron pairs, corresponding to the decay N 4 → νe + e − . As can be seen, for M 4 < m π the leading production mechanism is π ± decay. For masses in the region m π < M 4 < m K , K ± dominates and, in fact, the sensitivity contour reaches lower values of U e4 at the best point (in spite of the smaller number of kaons produced, when compared to the number of pions). The reason for this is that for M 4 < m π the heavy neutrino becomes very long-lived, leading to a reduced number of decays inside the detector and a consequent reduction in sensitivity. On the other hand, in the heavy mass region (M 4 > m K ) the heavy neutrino is predominantly produced from either D or D s meson decays (although there is a subdominant contribution from τ decays). While the D s is heavier (and therefore more difficult to produce) than D mesons, its decay to heavy neutrinos is mediated by the CKM element V cs instead of V cd . This compensates for the reduced meson production rate and, as a result, the sensitivity in this region is dominated by D s decays. Finally, the different slope as a function of M 4 for the D contribution is simply due to the fact that, unlike for the π, K and D s decays, the D meson production of N 4 is dominated by three-body decays instead of two-body (see Sec. 4). Fig. 6 shows the sensitivity contours in the M 4 − |U α4 | 2 plane at 90% CL, for different decay channels as indicated. The upper, middle and lower panels in the figure show the results assuming that the heavy neutrino mixes predominantly with the e, µ and τ sectors respectively.
Finally, Fig. 7 summarizes in blue the 90% CL expected sensitivities at the DUNE near detector to the heavy neutrino mixing |U α4 | 2 as a function of its mass, assuming a Dirac HNL. In the Majorana case, the increase in the number of events would translate into a slightly better sensitivity, although the results would be qualitatively very similar. In this last figure we combine the events from all the channels depicted in Fig. 6 under the same assumption of 20% signal efficiency and negligible background, following Ref. [84]. We again estimate the sensitivity following the Feldman and Cousins [85] prescription for a Poisson distribution under the hypothesis of no events being observed, which corresponds to the expected total number of signal events combining all channels leading to a visible final state in the detector being smaller than 2.44. Figure 6. Expected DUNE sensitivity (at 90% CL) to the mixing matrix elements |U α4 | 2 as a function of the heavy neutrino mass, for a total of 7.7 · 10 21 PoT collected. In each row, we assume that the HNL only couples to one of the charged leptons as indicated, while the other two mixings are set to zero. The different regions correspond to the results for different final states as indicated by the labels. Left panels correspond to signatures with charged leptons and missing energy, while middle (right) panels correspond to signatures with pseudoscalar (vector) mesons in the final state. In our analysis, we assume a negligible background level after cuts and a signal selection efficiency of 20%, see text for details. For comparison, the shaded gray areas indicate the parameter space disfavored by current experiments (at 90% CL). Relevant bounds on U e4 are obtained from results by the TRIUMF [86], PIENU [87], NA62 [88], T2K [84], CHARM [89], BEBC [90] and DELPHI [91] collaborations; for U µ4 , by PSI [92], PIENU [93], PS191 [94,95], 3 E949 [97], T2K [84], NuTeV [98] and DELPHI [91]; finally, U τ 4 is much harder to probe experimentally and here the only available constraints come from CHARM [99] and DELPHI [91] 4 . We find that DUNE is expected to improve over present constraints by several orders of magnitude in a large fraction of the parameter space and, in particular, for HNL masses between the K and D meson thresholds.
As a target region, we have also indicated in Fig. 7 the naive expectation for the mixing matrix elements from the Seesaw mechanism: |U α4 | 2 ∼ m i /M 4 , where m i stands for the SM neutrino masses. In particular, we set as the lower end of the band the minimum mass that at least one of the neutrinos must have to correctly reproduce the atmospheric mass splitting as measured in neutrino oscillations ∆m 2 atm = 0.05 eV. The upper line has been set using the latest bound of 1.1 eV from the KATRIN experiment [100]. We find that DUNE will be able to start exploring this interesting region, for HNL masses close to the K mass. Notice that this is only a generic expectation from the Seesaw mechanism: individual elements of the mixing matrix could either exceed or fall below these limits.
If sufficiently long lived, HNLs could decay during Big Bang Nucleosynthesis (BBN), altering the prediction for the primordial abundance of light elements. Thus, too small mixings are disfavoured, especially at low masses, as they would imply too long lived HNLs. BBN constraints exclude squared mixings smaller than ∼ 10 −5 for a neutrino with a mass of 100 MeV, while the bounds are much looser for larger masses, only disfavouring squared mixings below 10 −10 for an HNL with a mass of 1 GeV [101][102][103][104][105]. 5 Nevertheless, since these constraints rely on the cosmological history of the Universe, we choose not to display them together with direct laboratory tests in Fig. 7.
The number of HNL events depends on both their production rate and their decay probability inside the detector. At low masses the heavy neutrino production is dominated by pion decay, which is roughly proportional to |U α4 | 2 M 2 4 (see Eq. (4.1)). In this region, the most important HNL decay channel is N 4 → νe + e − , which is proportional to |U | 2 M 5 4 (see Eq. (5.8)). Thus, according to Eq. (6.3), the number of events should scale as |U α4 | 4 M 8 4 (an extra M 4 power arises due to the 1/γ factor, proportional to M 4 ). We have indeed verified that the slopes in the low mass regions of Fig. 7 fit well to |U α4 | 2 ∝ M −4 4 , as expected. We have also compared our results to similar studies in the literature [53,54,56], after the corresponding rescaling of the number of events accounting for the different detector volumes, PoT and efficiencies assumed, we find a rather good agreement between the four estimations of the DUNE sensitivity. We find the best overall agreement with Ref. [56]. The main difference is a slightly better sensitivity to the U τ 4 mixing in our results, in the small sensitivity peak we find around M 4 ∼ 1 GeV, corresponding to the closure of the τ → N 4 ρ − production channel. This peak also seems absent in the other references. Regarding Ref. [53], the main differences we find are at the peaks in sensitivity at the kinematic thresholds of the meson masses, where we find better sensitivity. We believe that these differences are due to the effect of the boost factor on the detector acceptance discussed in subsection 6.1, which becomes most relevant close to the kinematic thresholds, as shown in Fig. 4. We also find that the sensitivity to U µ4 for values of M 4 larger than the Kaon mass, is significantly smaller in Ref. [53] as compared to the other estimations, which find a similar behavior to that of U e4 , as expected from their similar branching ratios. Finally, we also find generally good agreement with Ref. [54]. The main differences are in the areas of parameter space were the HNL decays to ρ 0 and especially to π 0 are most relevant, since these decay modes were not included. The slope of the sensitivity curves is also slightly less steep than the |U α4 | 2 ∝ M −4 4 found in the other references.

Summary and conclusions
The addition of at least two nearly-sterile neutrinos (or HNLs) to the SM particle content is the simplest extension of the SM capable of reproducing the observed pattern of neutrino masses and mixing. The Majorana mass scale, unlike the masses of the other elementary particles, is not related to the electroweak scale and is a priori a free parameter of the model. The phenomenological consequences due to the existence of such heavy neutrinos would be very diverse depending on its value. In fact, while traditional type-I Seesaw

Neutral mesons
Charged mesons Pseudoscalars models set their Majorana masses at very high energies (experimentally inaccessible), lowerenergy versions (with heavy neutrinos at around the GeV scale) have recently drawn a lot of attention in the community since they are testable, do not worsen the hierarchy problem, and are able to reproduce the observed Baryon Asymmetry of the Universe. In such lowscale Seesaw models, the new singlets may form a pseudo-Dirac pair and lepton number is approximately preserved in the theory. The most promising avenues to look for MeV-to GeV-scale neutrinos are peak searches in meson decays, and searches for displaced vertices in fixed target experiments (produced when the neutrino travels a macroscopic distance before decaying back to SM particles). In both cases, an effective theory describing the interactions at low energies between mesons, neutrinos and charged leptons, obtained after the electroweak bosons have been integrated out, is the most suitable description. While most relevant vertices of the effective theory had been partially derived in previous literature, several inconsistencies remained. In this work, we have systematically derived all effective vertices involving mesons with masses of up to 2 GeV with significant branching ratios into HNLs. This allowed us to derive analytic expressions for the decay widths of the heavy neutrino into the different channels, and to clarify the inconsistencies found in previous literature (summarized in Sec. 5.4). For convenience, Tab. 7 summarizes the Feynman rules for the effective vertices involving charged leptons, mesons and neutrinos.
Our results have been made publicly available as FeynRules models [67] so that not only the total widths, but also fully differential event distributions, can be computed using Monte Carlo generators such as MadGraph5 [68]. This has been done separately for Dirac and Majorana HNLs. Moreover, note that, while the present work focuses on the lowenergy theory, our FeynRules implementation is more general and includes an option to replace all mesons with quarks, so they may also be used to study HNL phenomenology in collider searches at higher energies.
To illustrate the applicability of the effective theory and its FeynRules implementation, we have performed numerical simulations to obtain the expected heavy neutrino flux that would reach the DUNE near detector (ND), as well as the expected number of HNL decays inside the detector into several decay channels. The very high beam intensity, combined with the availability of a ND located at a distance L ∼ O(500) m, puts the DUNE experiment in an ideal position to search for the decay signals of HNLs produced from meson decays. We have shown how a proper treatment of the boost of the heavy neutrino, accounting for its mass, leads to an increased detector acceptance for the heavy neutrino flux when compared to the light neutrino case, see Figs. 3 and 4.
Finally, while the computation of the expected sensitivity at DUNE eventually needs a detailed detector simulation to address background rejection, it has been shown that, applying proper kinematic cuts to the particles observed in the final state, it is possible to reduce the background to a negligible level while keeping most of the signal events. Under this assumption, we have estimated in Sec. 6 the expected sensitivities to the model as a function of the heavy neutrino mass. We find that DUNE is expected to reach sensitivities comparable to or even better than those of fixed target experiments (see Fig. 7). We also find that DUNE will be expected to start exploring the region of parameter space where neutrino masses can be explained using a type-I Seesaw model, for HNL masses around the K mass scale.

A Generators of SU(3)
As outlined in Sec. 3, the normalization for the SU(3) generators has been chosen to satisfy the trace conditions Tr For convenience, we provide explicit expressions for the SU(3) generators below:

B Mixing matrices in a 3 + 1 scenario
It can be interesting to consider a case in which only one heavy neutrino is light enough or exhibits a sufficiently large mixing to play a role in the relevant phenomenology 6 . In this case, the model parameters will be four: the three leptonic Yukawa couplings Y ν,α and the heavy mass M , defined in Eqs. (2.1) or (2.3) (for the type-I and inverse Seesaw models, respectively). It is possible to write a 4 × 4 mixing matrix U in terms of these parameters, which rotates from the flavor basis to the mass one. The shape of such matrix will depend on whether neutrinos are either Majorana or Dirac fermions. For the Dirac case, the mixing matrix U relates the 4 left-handed neutrinos (ν L,e , ν L,µ , ν L,τ , N L ) to the 4 mass eigenstates (n 1 , n 2 , n 3 , N 4 ). In terms of the parameters mentioned above, the mixing matrix reads: where θ α ≡ Y ν,α v/ √ 2M , θ 2 ≡ |θ e | 2 + |θ µ | 2 + |θ τ | 2 and r ≡ √ 1 + θ 2 . In this case, only N 4 is massive, with a Dirac mass M 4 = rM . In the limit in which all the mixing parameters θ α are small, r ∼ 1 and the mass of the heavy neutrino will be approximately M . 6 The masses of the light neutrinos can be neglected for phenomenological purposes here.
On the other hand, in the Majorana case the flavor eigenstates are (ν L,e , ν L,µ , ν L,τ , N c R ). The mixing matrix now takes the form: with θ s ≡ θ e + θ µ + θ τ and ρ ≡ 1/ √ 1 + 4θ 2 . The i factors are chosen to obtain positive masses when diagonalizing the mass matrix. Now only two mass eigenstates, n 1 and n 2 , are massless, while n 3 and N 4 have Majorana masses of M 2 |1 ∓ ρ −1 | respectively. If all the mixing parameters θ α are small, then ρ ∼ 1 , so the mass of n 3 is negligible and that of N 4 is approximately equal to M .

C Determination of the vector meson decay constants
Unlike pseudoscalar mesons, vector meson resonances are wide and unstable under QCD. Thus, the determination of their decay constants is generally challenging, with more variability among different estimations in the literature. In order to bypass this issue, a possibility is to compute the width for a decay channel mediated by the electroweak interaction that has been precisely measured, comparing the result to the experimental values from Ref. [69]. This way the corresponding value of the decay constant can be directly extracted for each of the resonances under consideration, ensuring that the notation and normalization conventions used are consistent.

C.1 Neutral vector mesons
In this case, a good choice is the decay channel V → e + e − , which has been precisely measured and is dominated by photon exchange. Thus, we decompose the electromagnetic (EM) current j V EM,µ = i q eQ qq γ µ q as a linear combination of the meson currents, as we did for the Z current in Sec. 3: Comparing these results to the corresponding measurements [69,106], we find the values for the decay constants f V listed in Tab. 1.

C.2 Charged vector mesons
For the ρ ± mesons, we will use the ρ 0 constant, already determined, since the isospin breaking corrections should be negligible. However, for the K * ,± meson we must compute the decay width for an electroweak process, extracting the decay constant from there as we did for the neutral vector mesons. In this case, a good choice is the process τ − → K * ,− ν τ . The authors of Ref. [107] perform such calculation and report the value of the ratio between the ρ and K * decay constants: Therefore, using f ρ = 0.171 GeV 2 , we obtain f K * = 0.178 GeV 2 as listed in Tab. 1.

D Implementation of semileptonic form factors into FeynRules
The form factors involved in semileptonic meson decays include a dependence on the squared momentum transfer between the involved mesons, q 2 , which is not trivial to implement in a UFO model. For this reason, we have included two different implementation choices into our FeynRules models: a simpler option, which neglects the q 2 dependence and has been tuned to approximately reproduce the correct branching ratios for semileptonic decay channels; and a more sophisticated one, which includes the correct q 2 dependence as described in Sec. 3.3.1. As a first option, our FeynRules model files include by default constant form factors, evaluated at an average value of the (squared) momentum transfer, ⟨q 2 ⟩. This average value is determined by imposing that the correct total decay width is obtained, and depends mildly on the heavy neutrino mass M 4 and the charged lepton mass. We perform a cubic fit for the function ⟨q 2 ⟩(M 4 ) in the kinematically allowed range for the heavy neutrino mass. Furthermore, the coefficients of such a fit depend on the flavor of the involved charged lepton, in order to account for the correct kinematics. The accuracy of this approximation relies on the fact that the dependence of the form factors on the momentum transfer is very mild in the allowed kinematic range.
As a second option, together with the FeynRules model files we provide a Python script that, upon running, modifies the relevant files in the output UFO. This way, the correct energy dependence of the vertices, according to the linear and pole parametrizations of the form factors described in Sec. 3.3.1, can be implemented, allowing for precise event generation in MadGraph5.
In short, if the provided Python script is run after generating the UFO file with Feyn-Rules, the correct energy dependence of the form factors can be fully incorporated into MadGraph5. Otherwise, the constant form factors evaluated at ⟨q 2 ⟩ in the default Feyn-Rules model allow for a good approximation also for the other formats into which the model may be exported to.