Effective Majorana mass matrix from tau and pseudoscalar meson lepton number violating decays

An observation of any lepton number violating process will undoubtedly point towards the existence of new physics and indirectly to the clear Majorana nature of the exchanged fermion. In this work, we explore the potential of a minimal extension of the Standard Model via heavy sterile fermions with masses in the $[ 0.1 - 10]$ GeV range concerning an extensive array of"neutrinoless"meson and tau decay processes. We assume that the Majorana neutrinos are produced on-shell, and focus on three-body decays. We conduct an update on the bounds on the active-sterile mixing elements, $|U_{\ell_\alpha 4} U_{\ell_\beta 4}|$, taking into account the most recent experimental bounds (and constraints) and new theoretical inputs, as well as the effects of a finite detector, imposing that the heavy neutrino decay within the detector. This allows to establish up-to-date comprehensive constraints on the sterile fermion parameter space. Our results suggest that the branching fractions of several decays are close to current sensitivities (likely within reach of future facilities), some being already in conflict with current data (as is the case of $K^+ \to \ell_\alpha^+ \ell_\beta^+ \pi^-$, and $\tau^- \to \mu^+ \pi^- \pi^-$). We use these processes to extract constraints on all entries of an enlarged definition of a $3\times 3$"effective"Majorana neutrino mass matrix $m_{\nu}^{\alpha \beta}$.


Introduction
In parallel to the direct searches being currently carried at the LHC, the high-intensity frontier also offers a rich laboratory to look for (very) rare processes which indirectly manifest the presence of new physics. Individual lepton flavours and total lepton number are strictly conserved quantities in the Standard Model (SM). By themselves, neutrino oscillations constitute evidence for lepton flavour violation in the neutral lepton sector, and suggest the need to consider extensions of the SM capable of accounting for the necessarily massive neutrinos and for lepton mixing, as confirmed by experimental data [1][2][3][4][5][6][7][8].
The observation of a lepton number violation (LNV) process would be a clear signal of new physics, in particular of the existence of Majorana fermions. Although neutrinoless double beta decay (0ν2β) remains by excellence the observable associated with the existence of Majorana neutrinos, many other processes reflecting ∆L = 2 are being actively searched for. At colliders, there are several possible signatures of LNV and/or manifestations of the existence of Majorana fermions (see for instance [9][10][11][12][13][14] and references therein); in some cases these might even hint on the neutrino mass generation mechanism (if interpreted in the light of such theoretical frameworks).
Semileptonic meson decays are examples of transitions which offer the possibility of studying LNV, for both cases of three-or four-body final states, M 1 → ℓ ± ℓ ± M 2 and/or M 1 → ℓ ± ℓ ± M 2 M 3 , M i denoting mesons. Likewise, there are several semileptonic tau decays corresponding to ∆L = 2 transitions, τ ± → ℓ ∓ M 1 M 2 , which can occur in the presence of Majorana neutral fermions (such as neutrinos). However, whether such decays do have a non-negligible width strongly depends not only on the Majorana nature of the neutrinos, but on their properties. If the processes are only mediated by the three active light neutrinos, the corresponding widths will be strongly suppressed by the smallness of their masses (as they are proportional to m 2 ν ); should one consider extensions of the SM particle content by heavier neutral fermion states, which have non-vanishing mixings with the left-handed neutrinos of the SM, then the additional contribution to the widths will also be suppressed, either by the heavy propagator's mass, or by tiny mixings. The suppression can be nevertheless evaded if the mass of the heavy neutrino is such that it can be produced on-shell in the decay of the heavy meson (or of the tau lepton). Due to the so-called "resonant enhancement", the LNV decay widths can be strongly increased in this regime.
In recent years, the rôle of heavy sterile fermions (sufficiently light to be produced on-shell from the semileptonic LNV decays of mesons or taus) has been addressed, both for three-and four-body final states [9,11, and bounds were derived for the new propagator's degrees of freedom. Just as valuable information on the ee entry of the effective Majorana mass matrix can be inferred from LNV 0ν2β decays (see, e.g. [36]), bounds have been established for other entries of m αβ ν , using adapted expressions inspired from neutrinoless double beta decays [10,28,30,[37][38][39][40]. B-factories (such as BaBar and Belle), together with the advent of the LHC (in its run 1 and 2), have allowed to establish increasingly stronger limits on LNV semileptonic decays [41,42]. On the theoretical side, significant progress has been made on the non-perturbative computation of several quantities -such as decay constants -relevant for the LNV semileptonic transitions. Although these decays all imply the violation of individual lepton flavours (by two units), some can lead to the appearance of same sign but distinct (i.e., different flavour) charged leptons in the final state, M 1 → ℓ ± α ℓ ± β M 2 with α = β. At the high-intensity frontier, a number of experiments dedicated to search for a signal of charged lepton flavour violation (cLFV) -already collecting data, or due to operate in the near future -are expected to significantly improve the bounds on several processes. Such is the case of µ → eγ for which the present (future) bound is 4.2 × 10 −13 [43] (6 × 10 −14 [44]), of µ → eee with a present (future) bound of 1.0 × 10 −12 [45] (∼ 10 −16 [46]), and of µ − e conversion in nuclei with a present bound of 7 × 10 −13 for gold nuclei [47] and an expected sensitivity of ∼ 10 −17 [48][49][50] for Aluminium targets. Likewise, valuable information can also be obtained from cLFV meson decays, such as M → ℓ + α ℓ − β and M 1 → M 2 ℓ + α ℓ − β , currently searched for both at B-factories and at the LHC [41].
In view of the recent experimental and theoretical progress, and in preparation for the expected improvement in the experimental fronts, in this work we explore a variety of semileptonic meson and tau LNV decays induced by the presence of sterile Majorana fermions, whose mass is such as to allow for the resonant enhancement of the decay rates above referred to. Moreover, we propose a generalised definition of the effective Majorana mass matrix (encompassing the standard m ee ν associated with 0ν2β decays), and infer constraints on all its entries from the confrontation of a comprehensive set of LNV semileptonic decays to current experimental bounds.
In the present study, we only consider processes leading to three-body final states -in particular, when considering LNV meson decays, we only address processes involving a single meson in the final state; furthermore, we focus on pseudoscalar meson decays. This implies that we will not address channels such as the decays of vector mesons, nor the LNV four-body decays of B mesons, which have recently been identified as promising experimental probes, in particular at LHCb (see, e.g. [12,21,[32][33][34][35]). Despite their experimental appeal, these decays are still subject to important theoretical uncertainties, such as poor control of the hadronic matrix elements (when vector mesons are involved), or interactions between final state hadrons, among others; we will not address them in the present study.
Concerning the underlying theoretical framework, and as a first step, we do not consider a specific mechanism of neutrino mass generation, but rather a minimal simple extension of the SM via one Majorana sterile fermion, focusing on the [0.1 − 10] GeV mass range. Any scenario of new physics involving sterile fermions -even such minimal constructions -must comply with stringent phenomenological and observational constraints arising from bounds derived from different experimental frontiers (high-intensity, high-energy experiments, neutrino data, as well as from cosmology). All available constraints will be taken into account in our analysis.
We conduct a systematic analysis of lepton number violating semileptonic meson decays, also including lepton number violating tau-lepton decays. This allows to build upon and update earlier analysis, revising the bounds derived in [10,11,28,30,[37][38][39][40][51][52][53], and obtaining an up-to-date overview of the sterile fermion parameter space (mass and combinations of the active-sterile mixing angles, U ℓα4 U ℓ β 4 ). We explore the allowed parameter space in order to identify regimes for which the LNV branching ratios might be within experimental reach, for instance of NA62 (for the light mesons), BES-III for charmed mesons and LHCb (and in the future, Belle II) for B mesons and tau leptons. (Our study identified decays which already are in conflict with current bounds, as is the case of K + → ℓ + α ℓ + β π − , and τ − → µ + π − π − .) We also impose that the on-shell heavy fermion be sufficiently short-lived as to decay within a finite size detector (we will work for a "benchmark" value of 10 metres).
As mentioned above, we further translate all the collected bounds into limits for the flavour conserving and flavour violating entries of the 3 × 3 Majorana effective mass matrix, |m αβ ν |. Other than (strongly) improving existing bounds, we propose bounds for the |m τ τ ν | entry. This work is organised as follows: in Section 2 we update the constraints on minimal SM extensions by a sterile neutral fermion arising from LNV (and cLFV) semileptonic three-body decays, reviewing their experimental status, and we revisit the theoretical estimation of the corresponding decay widths in the "resonant enhancement" regime; once all phenomenological and observational constraints have been applied, we update the sterile fermion parameter space. We further address the reconstruction of the 3 × 3 "effective" Majorana mass matrix. Section 3 collects our results (and predictions): LNV branching fractions for the semileptonic decays of mesons (B, D, K) and τ leptons, as well as their impact for the effective Majorana mass entries. Further discussion and remarks are given in Section 4. Other than the detailed description of the minimal SM extension we consider, Appendix A summarises the relevant constraints applied in our work. Appendix B details the computation of the LNV meson and tau lepton decay widths, while in Appendix C we describe the calculation of the decay width of the exchanged Majorana state. Finally, the expressions for the cLFV meson decay widths are collected in Appendix D.
2 An additional sterile fermion: impact for LNV and cLFV Extensions of the SM via the addition of Majorana sterile neutral fermions, with non-negligible mixings with the light, active neutrinos, open the door to numerous processes which violate total lepton number, charged lepton flavour (or both), among other phenomena. While many complete models of neutrino mass generation do include such states (for example right-handed neutrinos in several seesaw realisations), considering a minimal scenario in which a single Majorana sterile neutrino is added to the SM field content 1 , without any assumption on the mechanism of neutrino mass generation, proves to be a useful first step in evaluating the potential contribution of these states to a wide array of observables, including cLFV and LNV decays.
The aim of this section is to discuss the updated constraints on minimal SM extensions by sterile fermions arising from LNV (and cLFV) semileptonic three-body decays. After a summary of the corresponding experimental status, we discuss the theoretical computation of the decay widths, and update the constraints on SM extensions with sterile Majorana fermions from the non-observation of LNV decays. We further address the reconstruction of the effective 3 × 3 Majorana mass matrix. In order to do so, we thus rely on a simple "toy model" which is built under the single hypothesis that interaction and physical neutrino eigenstates are related via a 4 × 4 unitary mixing matrix, U ij (see Appendix A). Other than the masses of the three light (mostly active) neutrinos, and their mixing parameters, the simple "3+1 model" is parametrised by the heavier (mostly sterile) neutrino mass m 4 , three active-sterile mixing angles, as well as three new CP violating phases (two Dirac and one Majorana).
In our discussion, and unless otherwise stated, the sterile state is assumed to be produced (on-shell) from the semileptonic decays of either tau leptons or mesons.

LNV and cLFV decays: experimental status
As mentioned in the Introduction, searches for LNV and/or cLFV decays are (and have been) an important goal of numerous experimental collaborations.
Especially in a context in which LNV arises from the presence of Majorana fermions, neutrinoless double beta decay is one of the most promising observables: KamLAND-ZEN [54,55], GERDA [56] and EXO-200 [57][58][59], have all set strong bounds on the m ee ν effective mass, to which the amplitude of 0ν2β process is proportional (see detailed expressions in Section 2.4). The most recent and strongest constraint has been obtained by the KamLAND-ZEN collaboration [55] m ee ν 0.165 eV (90% C.L.) ; (1) concerning future prospects, we summarise in Table 1 the sensitivity of ongoing and planned 0ν2β dedicated experiments. In this work, we take a representative benchmark for the future sensitivity |m ee ν | 0.01 eV (for details concerning the theoretical uncertainties, see for instance [60,61]).

Experiment
Ref.
Semileptonic meson decays offer numerous channels to look for LNV; final states comprising two same-sign leptons and one or two mesons, M 1 → ℓ ± ℓ ± M 2 and M 1 → ℓ ± ℓ ± M 2 M 3 , have been experimentally searched for in recent years. In Table 2 we collect the bounds for the case of LNV three-body final states (including same and different charged leptons in the final states), further referring to the m αβ ν effective mass entries thus constrained. The current bounds listed are expected to be improved by one or two orders of magnitude with future experiments, such as LHCb in its run 2 [75] as well as NA62 [76]; for the channel B + → µ + µ + π − , the future experiments Belle II 2 and FCC-ee are expected to improve the current bounds (for |U µ4 | 2 ) by about two orders of magnitude [31]. Although not included in Table 2, the decay mode B − c → τ − τ − π + can also be explored aiming at constraining the matrix element m τ τ ν [26] (and other elements [17,27] as well).
Important constraints, albeit indirect, on the effective mass entries m αβ ν can also be inferred from lepton number conserving (but cLFV) semileptonic meson decays, as the latter constrain individual entries (or combinations) of the lepton mixing matrix, U ℓαi U ℓ β i . While cLFV radiative and three-body decays (such as ℓ α → ℓ β γ, ℓ α → 3ℓ β , or neutrinoless conversion in Nuclei) in general lead to severe constraints for sterile masses typically above the tenths of GeV, in the mass regimes in which the heavy neutrino is produced on-shell from the above discussed tau and/or meson decays, constraints arising from cLFV meson decays can be important, and should thus be carefully evaluated and taken into account. In Tables 3 and 4 we summarise some of these bounds, which are of relevance to our study.
Due to its higher mass, the tau lepton can also decay semileptonically and, in the presence of new physics sources several channels might signal ∆L = 2 transitions such as τ ± → ℓ ∓ M 1 M 2 (which is necessarily cLFV). The current bounds regarding τ decays are summarised in Table 5, and are used to infer constraints on the m ℓτ ν effective mass entries.

Meson and tau lepton decay widths
We now proceed to discuss and highlight relevant points leading to the computation (theoretical derivation) of the LNV meson and tau semileptonic decay widths. These are schematically depicted in Fig. 1 for the case of a semileptonic LNV meson decay.

Theoretical estimation
As already mentioned, leading to the computation of the LNV semileptonic decays, we have made several assumptions, which we proceed to discuss.
• We consider semileptonic decay modes leading to three-body final states; moreover, we only consider the decays of pseudoscalar mesons and do not address vector meson decays, as their (non-perturbative) decay constants are plagued by larger theoretical uncertainties, and the resonances (and excitations) are not well determined; • The only source of lepton number violation (and lepton flavour violation) at the origin of the distinct decays above mentioned stems from the presence of (heavy) Majorana neutrinos; • In order to avoid excessive suppression due to the propagation of a virtual heavy state, we focus on the resonant regimes in which, depending on the decaying particle and final state considered, sterile neutrinos with a mass within the kinematically allowed interval are produced on-shell, and propagate before decaying; • Via their interaction with the light (mostly active) neutrinos, the extra sterile fermions lead to a modification of the leptonic neutral and charged currents, as given by Eqs. (11), (13) of Appendix A.1. In turn, this means that many other observables are also affected, hence leading to extensive experimental and observational constraints on the sterile neutrino degrees of freedom; the latter are summarised in Appendix A.3.
The full expressions for the LNV decay widths of mesons and tau leptons 3 can be found in Appendix B; they can be however cast in the following compact form, in which V M denotes the CKM matrix element relevant for the quark content of meson M , and f M is the meson's decay constant. Likewise, U ℓi is the leptonic mixing matrix element for the W ℓν interaction. The function Φ encodes the integrals over the corresponding phase space, and depends on the masses and momenta involved in the transitions. As mentioned before, the decay widths are computed under the assumption of an on-shell neutrino; moreover, and as detailed in Appendix B, one works in the narrow-width approximation, which is valid and appropriate in the regimes here discussed. In recent years, important progress has been made regarding meson data, particularly in what concerns the computation of the decay constants which are crucial in the estimation of Γ LNV (see Eq. (2)). In Tables 6 -10 we summarise the most recent data (the reported values are either experimentally determined, or non-perturbatively computed) for the mesons relevant to our analysis.  Table 6: Decay constants and masses for pseudoscalar mesons: π ± , K ± , D ± , D ± s , B ± , as well as η, η ′ and η c .   [86].
Although in our study we do not address the decays of vector mesons, the corresponding decay constants and masses are collected in Tables 8, 9 and 10, as such inputs are relevant regarding the computation of the (on-shell) heavy neutrino decay width.   [93].  Table 9: Decay constants and masses for ρ, K * and D * (s) vector meson states (the values for the decay constants of K * 0 ,K * 0 (D * 0 ,D * 0 ) are the same as those of K * ± (D * ± ) since they are related via isospin symmetry).

On-shell sterile neutrinos and finite detector effects
As previously mentioned, LNV semileptonic meson and tau decays can be strongly enhanced when the exchanged Majorana fermion is on-shell. This is equivalent to having the production of a real state, which propagates and has a well-defined width. Depending on its mixing with the leptonic doublets, a sterile state with a mass in the range relevant for our study (i.e. 100 MeV m νs 10 GeV) can have several possible decay channels. In addition to those which are directly related with the LNV final states, ν s → ℓM , with M denoting a meson, it can also have (simple) twobody decays, ν s → ν i ℓ, i = 1, 2, 3. In the present study, we carry a full computation of the heavy neutrino decay width, including all possible channels which are open for a given mass range. We collect the relevant expressions in Appendix C.
It is important to stress that the observability of the different processes referred to in the previous section is not only related to the expected number of events, but also to whether or not the final states can be indeed observed -if sufficiently long-lived, the LNV decays of the heavy neutrino can occur outside the detector. Thus, in our analysis we will impose that the distance travelled prior to the decay does not exceed the typical size of a detector. The actual distance travelled by the neutrino depends on its lifetime and on its velocity β s , L flight νs = γ s β s τ s c. For relatively heavy neutrinos (with a mass typically above 500 MeV), τ s can be sufficiently small to ensure that the length of flight does not exceed a few metres (or even less). On the other hand, for lighter states, the mixings tend to be very small (due to the important experimental constraints, see Appendix. A.3), and such states might decay well outside the detector. In the latter case, and independently of the actual decay widths, the final states would not be observed. In Fig. 2 we display the heavy neutrino length of flight, obtained for the maximal allowed mixings to the active states. This corresponds to a lower bound on L flight νs , as the lifetime will be longer for smaller mixings. For instance, a heavy neutrino with a mass around 150 MeV produced with an energy of ≃ 200 MeV, already has γ s β s ∼ 0.88, and will certainly escape detection.
In order to have realistic results concerning the constraints arising from the LNV bounds, the effect of a real (finite) detector volume must then be taken into account upon computation of the LNV branching ratios: the probability that the sterile neutrino length of flight is smaller than the detector size, L flight νs L det. , is given by where p νs , m νs and Γ νs are the momentum, mass and decay width of the heavy (mostly) sterile neutrino, respectively. The momentum p νs depends on the energy of the mother particle (M 1 or τ ) produced by experiments; as done in Ref. [15], here we simply consider that the mother particle is produced at rest or with a comparatively small momentum. Should that be the case, p νs is then given by where m A = m M 1 and m B = m ℓα , m ℓ β respectively for the LNV meson decay, or m A = m τ and m B = m M 1 , m M 2 for the LNV τ decay; m νs is the mass of the sterile state ν s and the kinematical function λ(m 2 A , m 2 νs , m 2 B ) is defined in Eq. (44). The theoretically computed branching fraction of each LNV decay channel is then corrected by the within-detector decay probability (cf. Eq. (3)).
In what follows -and as already stated in the beginning of this section -we will consider a simple extension of the SM by one sterile Majorana fermion, which mixes with the light (mostly active) neutrinos. We will thus subsequently denote the corresponding parameters as m νs = m 4 , p νs = p 4 , Γ νs = Γ 4 , τ s = τ 4 and L flight νs = L flight ν 4 .

Updated constraints from LNV decays
Due to the significant contributions of a wide array of experiments, the sterile neutrino parameter space (as spanned by its degrees of freedom, m 4 and U ℓα4 , ℓ α = e, µ, τ ) is already subject to extensive constraints. Before revisiting those arising from LNV, we summarise in Fig. 3 the current status of the sterile neutrino parameter space (the bounds here imposed are enumerated and described in Appendix A.3). The different panels displaying the bounds on the active-sterile mixing angles (as a function of the mass of the heavy, mostly sterile neutrino) arise from direct searches, cLFV bounds, invisible Z decays, W → ℓν decays, among others [96][97][98][99][100][101][102][103]. Although not included in these plots, it is important to emphasise that leptonic observables such as R τ , R π and R K (or their corresponding deviations from the SM predictions, ∆r K , ∆r π and ∆r τ ) are extremely constraining, and will be taken into account in the subsequent numerical analyses.
Leading to the exclusion regions of Fig. 3, the limits on a given mixing element U ℓ4 are obtained by setting the other mixings with the sterile state to zero (for instance, the constraints on U µ4 are inferred assuming U e4 = U τ 4 = 0); for the exclusions on combined mixing elements U ℓ4 U ℓ ′ 4 , with ℓ = ℓ ′ , one assumes the latter to be equal, and the remaining one is set to zero (for instance, leading to the constraints on U e4 U µ4 , one has U e4 = U µ4 and U τ 4 = 0); finally, for the final panel regarding bounds on the triple product of couplings 4 (|U e4 U µ4 U τ 4 |), all entries are taken to be equal. The dot-dashed line on the first panel (|U e4 | 2 vs. m 4 ), corresponds to neutrinoless double beta decay searches, and is obtained by requiring that, by itself, the contribution of the sterile neutrino does not exceed 0.165 eV (we assume that no cancellations occur between the heavy and the light neutrino contributions); this should be interpreted as a representative bound, since it can be strengthened (relaxed) in the case of constructive (destructive) interference with the active neutrino contributions. The discontinuity in the exclusion for |U τ 4 | 2 vs. m 4 (for the interval m 4 ∈ [0.3, 0.5] GeV) arises from the lack of experimental bounds applicable to that particular mass regime. For completeness, the sterile mass regime constrained by the results of Fig. 3 extends to masses heavier than the mass interval characteristic of the meson and tau LNV decays here addressed (below 10 GeV); this is done in order to explicitly obtain a comprehensive picture of all available bounds -including those from cLFV decays (relevant for the combined mixing elements), which only become significantly constraining for sterile mass close to the electroweak scale (as manifest in the |U e4 U µ4 | vs. m 4 panel, for example). Bounds from invisible (leptonic) decays of Z (W ) bosons are also only relevant for a heavy sterile mass regime, and are typically superseded by the constraints arising from direct searches -and in the case of |U e4 | 2 , from neutrinoless double beta decay. Finally, the unitarity bound arising from non-standard neutrino interactions with matter [104] or from non-standard oscillation schemes, constrains the deviation from unitarity of the left-handed lepton mixing matrix (Ũ PMNS ), which leads to the exclusion of the large mixing regimes (for all the sterile mass regimes considered). In deriving these bounds we used the recent results [105] (see also [106]).
The different panels of Fig. 3 thus summarise the most recent and up-to-date available constraints on the parameter space of the SM extended by one additional massive Majorana fermion with a mass between 10 MeV and 100 GeV.
The bounds and exclusion regions on the sterile neutrino degrees of freedom, as identified in the different panels of Fig. 3, must now be combined with the bounds arising from LNV decays; together, they will have a strong impact on the maximal currently allowed values of the LNV semileptonic tau and meson decays, and most importantly, on the values of the entries of the 3 × 3 Majorana effective mass matrix.
We thus begin by discussing the updated constraints on the relevant combination of leptonic mixing matrix elements, stemming from the most recent experimental bounds. We recall that resonant production of sterile neutrinos (which is crucial to enhance the LNV meson and tau decays, implies that the sterile state is on-shell. Should the decay width of the sterile neutrino be sufficiently small, the decay may occur after . Solid surfaces denote excluded regimes due to violation of at least one experimental or observational bound. a length of flight larger than a realistic detector size, thus rendering the LNV decay processes invisible. Whenever relevant, the derived bounds take into account the requirement that the heavy neutrino decays within a finite detector (i.e. imposing that its length of flight does not exceed a nominal value L flight ν 4 = 10 m), as described in Section 2.2.2. The final results are then compared to the experimental limits listed in Tables 2 and 5. Working under the same assumptions as those leading to Fig. 3, we thus display in Figs. 4 and 5 the bounds on active-sterile mixing angles, as a function of the (mostly) sterile heavy neutrino mass, arising from LNV meson and tau decays. When present, dashed lines denote the bounds derived taking into account the requirement of having the heavy neutrino decaying within 10 m (see discussion above). These constraints will be subsequently used in our numerical analysis. As expected, given the current experimental bounds, the most stringent LNV constraints on the active-sterile mixing angles arise from semileptonic kaon decays (K + → ℓ + α ℓ + β π − ), leading to constraints on (combinations of) mixings of O(10 −9 ). Even when corrected to account for a finite detector size (within-detector decay), semileptonic kaon decays -especially leading to a final state containing at least one electron -are still the most stringent ones.
In recent years, similar analyses have led to the derivation of increasingly stronger bounds on the sterile fermion parameter space. Although very recent works already include updated experimental bounds in their results, our study considers the most recent data for all the LNV Lepton number violating semileptonic tau decays have a strong impact on constraining combinations of mixings involving the τ lepton, as can be inferred from a direct comparison of Fig. 3 with the corresponding panels of Fig. 5. It is interesting to notice that (as will be discussed in greater detail in Section 3.1.4) LNV tau decays lead to constraints on µ − τ mixing elements that can be as stringent as |U µ4 U τ 4 | 10 −5 for sterile mass close to 1 GeV -for conservative limits on the size of the detector. Such bounds already appear to supersede those displayed in the corresponding panel of Fig. 3.
The results collected in the above plots reflect two distinct regimes for the sterile neutrino lifetime. We have also investigated more extreme limits (for example L flight ν 4 1 m, as will be the case at Belle II), and our findings suggest a loss in sensitivity amounting to a factor 2 to 3 in the different combinations |U ℓα4 U ℓ β 4 |.
In order to compare our results with previous analyses carried in the literature, it is convenient to recast the current experimental bounds on LNV decays using a most minimal assumption for the active-sterile mixings, in particular that of degenerate mixing angles. Working under the hypothesis of |U e4 | = |U µ4 | = |U τ 4 |, the bounds on active-sterile mixing angles, as a function of the (mostly) sterile heavy neutrino mass, are displayed in Figs. 6 and 7.
We have thus compared our results with those obtained by other groups; as an example let us mention that our findings qualitatively agree with those of [11,15,35]. The general shape of the derived curves is in good agreement 5 , but our results lead to stronger exclusions, as they were obtained for the most recent experimental data. Moreover, we also emphasise that we do compute the full contributions to the on-shell heavy neutrino decay width.
Four-body semileptonic LNV meson decays 6 have also been explored in recent years as complementary means to constrain the active-sterile neutrino mixings in the kinematically allowed heavy neutrino mass intervals. In [21], the decays were first used to derive bounds on the mass of the on-shell Majorana neutrino. The LNV decays B → D ( * ) µ ± µ ± π ∓ (as well as B → µ ± µ ± π ∓ ) were addressed in [29,33]; the future sensitivity of LHCb might allow to improve upon the constraints on U µ4 for the relevant mass intervals. Very recent studies [34] further suggested that the 4-body decay modes of B 0 s → P − µ + µ + π − (P = K, D s ) could lead to a small amelioration of the limits inferred by LHCb from the LNV three-body decay B − → µ − µ − π + . The comparison of our results with the bounds derived for LNV 4-body B 0 decays in [35] -which appeared during the final stages of this work -suggest that while qualitatively akin in their constraining power, certain 3 body processes, such as K → eµπ or D s → µµπ, are in fact more stringent. Figure 6: Updated constraints on the relevant combination of leptonic mixing matrix elements (|U ℓα4 U ℓ β 4 |) arising from LNV pseudoscalar meson decays, as a function of the heavy sterile neutrino mass (GeV). Leading to the above plots, all active-sterile mixing elements were taken to be equal (i.e., |U e4 | = |U µ4 | = |U τ 4 |). Line and colour code as in Fig. 4.
Before proceeding to the following section we nevertheless re-emphasise that in general we will not work under the assumption of degenerate mixings. We also recall that we only consider three-body LNV final states. The general amplitude for a LNV process mediated by Majorana (active or sterile) neutrinos typically depends on the following elements: where q 2 is the momentum of the neutrinos mediating the process and Γ j is their corresponding decay width. This general expression can be simplified in limiting cases (depending on the properties of the exchanged Majorana particle): for the active neutrinos Γ i ≃ 0, m 2 i /q 2 ≪ 1 and one recovers the usual definition of the 0ν2β decay effective mass m ee ν ; if q 2 < 0 the (virtual) neutrinos are not on-shell, and Γ i can be neglected for the heavy states as well. Notice that in the presence of heavy (almost sterile) Majorana neutrinos, the 0ν2β effective mass is correctly parametrised by m ee ν −(125 MeV) 2 [107], in which q 2 = −(125 MeV) 2 is an average of the virtual momenta in different decaying nuclei.
Under the hypothesis that Majorana neutrinos are at the origin of all lepton number violating processes, the LNV meson semileptonic decays M 1 → ℓ ± α ℓ ± β M 2 allow to infer the corresponding contributions to the (α, β) entry of the "effective neutrino mass matrix", m αβ ν (q 2 ). The amplitude of the LNV decay processes under study (see Appendix B) includes the following terms in which q 2 > 0 denotes the momentum transfer, and Γ 4 the width of the additional sterile fermion with mass m 4 . For the specific LNV decay M 1 → ℓ ± α ℓ ± β M 2 , the transfer momentum is of order . The above equation thus becomes allowing to cast a generic definition of the neutrino effective mass matrix as (following the approach of [10]) The combination of lepton matrix elements entering in the computation of the effective mass is heavily constrained -not only from the LNV bounds, but from a plethora of other processes which, as mentioned above, includes cLFV decays, direct and indirect searches, as well as cosmology (see Section 2.3 and Appendix A.3). Although Eq. (7) does offer a means to infer constraints on the effective mass matrix, it is important to stress that the latter depend on the momentum transfer (q 2 , implicitly entering p 2 12 ), and thus no absolute bounds can be obtained 7 .

Results in minimal frameworks
We proceed to report the results obtained for the distinct types of LNV decays we have considered, emphasising on how the latter can be further constrained by considering the improved sensitivity to other modes also violating lepton flavour, and finally discussing their implication towards the constraining of the effective mass matrix. As mentioned in Section 2, we work in the framework of a minimal toy model (a "3+1" SM extension), described in Appendix A.2. Here we just recall that the simple "3+1 model" is parametrised by the heavier (mostly sterile) neutrino mass m 4 , three active-sterile mixing angles as well as three new CP violating phases (two Dirac and one Majorana). In the numerical analyses (which will in general assume a normal ordering for the light neutrino spectra) we will consider a range for the mass of the additional heavy state so that it can be produced as an on-shell state from the tau and meson decays, and we sample the active-sterile mixing angles from the interval [0, 2π] (likewise for the various CP violating phases).

LNV semileptonic decays
We begin our presentation of the numerical results by investigating the contributions of the onshell mostly sterile state concerning several lepton number violating processes. In particular, we investigate all kinematically viable 3-body decay modes of charged B, B c , D s , D and K mesons, as well as the 3-body LNV decay modes of the tau lepton. We consider both the semileptonic decays into pseudoscalar and vector meson final states.
Leading to the plots presented in this section -with the exception of those of Fig. 12 -in order to optimally enhance the LNV rates, one considers the maximally allowed active-sterile mixing regimes for the lepton flavours involved in a given decay, as obtained from the study whose results are summarised in Fig. 3. This allows to maximise the couplings responsible for the process, while setting to zero the mixings involving the flavour(s) absent from the decay (mother particle or products), which in turn minimises the neutrino decay width. In some mass regimes this mixing 7 Notice that in general this is indeed the case. Although in the usual (0ν2β) definition m ee ν is apparently independent of q 2 , this is a consequence of only having contributions of active (light) neutrinos (for which m 2 i ≪ q 2 ). In the presence of heavy (virtual) neutrinos, the momentum dependence is restored (although in a milder way with respect to the case of resonant enhancements considered in this work). pattern leads to a neutrino lifetime τ 4 such that τ 4 c > 10 m, rendering the process difficult to observe in realistic detectors. In these cases we provide an additional prediction for the expected branching fractions, obtained by considering non-zero mixings also with the flavours not involved in the process, in order to comply with the condition τ 4 c = 10 m. The presentation of our results (see e.g. Fig. 8 and similar ones in the present section) includes both predictions for the maximal allowed LNV branching ratios, encoding via a full line the more general result and using a dashed one whenever the condition τ 4 c < 10 m is taken into account. In the mass regimes where the two lines merge (or if no dashed line is present) the above condition is automatically fulfilled.
An important constraint in the mass range m 4 ∈ [0.14 − 0.49] GeV stems from the strong bound on cLFV kaon decays, BR(K + → π + µ + e − ) < 1.3 × 10 −11 [41], which can be mediated by on-shell neutrinos in this mass regime; by limiting the product U e4 U µ4 and/or the sterile neutrino width, this bound indirectly constrains the widths of the LNV channels involving both an electron and a muon as final states. When our default mixing pattern violates the experimental bound on K + → π + µ + e − we consider non-vanishing U τ 4 , thus allowing to comply with the latter bound. We have verified that the existing bounds on the remaining cLFV meson decay channels are not competitive with the constraints presented in Fig. 3.

B meson decays
The allowed branching fractions for the decays of the pseudoscalar B meson into pairs of leptons of same or different flavour and a meson are summarised in the panels of Figs. 8 and 9, for the case of pseudoscalar and vector meson final states. With the exception of the final states including a pair of tau leptons (in which only light mesons as K and π can be produced), all other decay modes include pions, kaons, as well as D and D s pseudoscalar mesons. For the case of vector mesons, ρ, K * , D * and D * s are kinematically allowed (with the exception of the τ τ mode, for which only ρ and K * are kinematically allowed). Whenever the regime would allow for long-lived Majorana mediators, we further recompute the maximal branching fractions imposing that the on-shell neutrino lifetime should comply with cτ 4 ≤ 10 m (dashed lines).
A common feature 8 to several of the curves displayed in Figs. 8 and 9 (as well as to others presented in this section) is the appearance of two regimes in the decay widths: while in some cases this is perceived as a double plateau, or a "bump" within a curve, in others it is evidenced as a clear discontinuity for given intervals of m 4 . To understand this, recall that the total decay amplitude stems from two contributions -the one directly displayed in Fig. 1, and the one corresponding to the exchange of the two lepton flavours as final states (either arising from the production or from the decay of the heavy neutrino). Under the assumption that one cannot distinguish the lepton's production vertex, the condition of having an on-shell heavy neutrino mediating, for instance, the process M 1 → ℓ α ℓ β M 2 (i.e. m M 1 > m 4 + m ℓα and m 4 > m ℓ β + m M 2 ), must be verified for both channels. Although in general the regions fulfilling the above conditions tend to overlap under the exchange α ↔ β (nevertheless giving rise to the double-plateau behaviour), this is not always the case, especially when m ℓα ≫ m ℓ β . An example can be found in B − → e + τ + D − s decay, where the on-shell mass regimes, [1.97, 3.50] GeV and [3.75, 5.28] GeV do not overlap. The contribution to the amplitude of the two regimes can also be very different -one regime might dominate over the other by as much as one order of magnitude, giving rise to "kinks" at the sterile neutrino mass values where the dominant contribution becomes kinematically forbidden. 8 Notice that for sterile mass regimes close to 0.1 GeV, the contour of the line associated with final states comprising a pion strongly reflects the excluded regions due to laboratory constraints (direct searches), as seen in   Table 2). From top to bottom, left to right, the final state lepton content is e + e + , µ + µ + , τ + τ + , e + µ + , e + τ + and µ + τ + .
Whether or not these processes are likely to be experimentally observed in the near futurefor example at the LHC run 2, be it at LHCb, CMS, or ATLAS -calls upon a more dedicated discussion, taking into account not only the ∆L = 2 decay widths, but the B-meson production prospects, and detector capabilities for such final states. A naïve approximation suggests the following estimated number of events, for example regarding B + → ℓ + α ℓ + β M − : in which, and for a given experiment, L int denotes the integrated luminosity, σ prod is the B meson production cross section, and D = ε D × ρ N corresponds to a factor which parametrises the overall detection efficiency for a given process (ε D ), times the acceptance factor of the detector (ρ N ). For the LHC operating at √ s = 13 TeV, the B meson production cross section has been recently reported by LHCb [108] to be σ prod LHCb (pp → B + + X) = (86.6 ± 0.5 ± 5.4 ± 3.4) µb, for L int LHCb = 0.3 fb −1 . Assuming that D LHCb should not be far from the percent level (i.e., ∼ 1%), one finds for the expected number of events Should the integrated luminosity augment to 10 fb −1 (or possibly to 50 fb −1 at the end of LHC run 3), then ∆L = 2 processes with branching fractions in the ballpark of O(10 −8 − 10 −10 ) could lead to a non-negligible number of events 9 . In particular, and given that the current bounds on B + → µ + µ + π − are already O(10 −9 ) (cf. Table 2), a future improvement of the experimental sensitivity 10 should be sensitive to ∆L = 2 transitions as induced by the SM minimally extended by a sterile fermion: as seen from Fig. 8 maximal branching fractions for the latter decay lie around O(10 −10 ), for sterile mass between 2 and 3 GeV. Naturally, this is also subject to having a conservative estimation of the detector size: assuming a smaller effective detector length might lead to different constraints on the sterile fermion parameter space, and thus modify the above discussion. For instance, for the case of a 1 m detector (as will be the case at Belle II) -which was already qualitatively discussed in what concerns the bounds on the combinations of the active-sterile mixing elements (see Section 2.3) -we have verified that the LNV branching rates would be reduced by between one and two orders of magnitude, depending on the actual regime of m 4 . Heavy mass regimes, as visible from the panels of Fig. 8, would not be sensitive to finite (realistic) detector effects.
In all analogy, the panels of Fig. 10 display the maximal (allowed) values of the LNV decays B + c → ℓ + α ℓ + β M − (final state pseudoscalar meson). With the exception of final states including at least one tau lepton (for which the final B − is kinematically inaccessible) all other decay modes include pions, kaons, D and D s , as well as B mesons. In the absence of experimental bounds (as no data is available concerning these very rare B c decays), one cannot comment upon the experimental impact of the LNV decays; however, our analysis suggests that the maximal associated branching fractions are larger than those obtained for B + mesons, by as much as two orders of magnitude. For example, for a sterile state with mass Figure 11 summarises the same information, but this time for final states including one vector meson.
Before addressing other decays, and for completeness, we consider the question of whether or not current cLFV bounds directly constrain the maximal allowed values of the LNV decay rates (leading to final states with leptons of same or different flavour). Recall that the most important contributions to cLFV observables typically emerge for sterile states with masses well above the 9 In [34], and working for similar operating benchmarks (luminosity and detector performance), it was found that a significant sensitivity could also be expected for LNV 4-body decays, in particular for B 0 s → K − π − µ + µ + and B 0 s → D − s π − µ + µ + . 10 In [33], the number of B mesons expected to be produced at the LHCb upgrade and at Belle II were reported to be 4.8 × 10 12 and 5 × 10 10 , respectively.  Table 2).
GeV -in particular for states heavier than the electroweak scale [109-111, 113, 114]; hence, and as can be verified from the panels of Fig. 12 (we chose semileptonic B decays into a pseudoscalar meson final state to illustrate this point), despite the expected correlation between the similar "flavour-content" transitions, cLFV bounds (current as well as the expected future sensitivities), do not directly constrain the LNV modes. We notice that leading to the results of Fig. 12, a more comprehensive survey of the parameter space has been conducted: the active-sterile mixings (both angles and phases) are randomly sampled (without any underlying hypothesis) from the It is relevant to notice that following our study of the impact of the sterile fermion on LNV semileptonic decays, a small subset of the points (denoted by green crosses in Fig. 12) is excluded due to being associated with excessive LNV decays, already in conflict with current bounds. This reinforces the expected experimental prospects of these observables.

D meson decays
We continue our discussion of LNV semileptonic meson decays by investigating the prospects for D s and D decays. For the latter decays, only pions and kaons (ρ and K * ) are candidates to the kinematically allowed final state pseudoscalar (vector) mesons. The results are summarised in the panels of Fig. 13. It is worth noticing that the maximal allowed rates for D + s → µ + µ + π − are very close to the current experimental bound; in the near future, this observable can thus play an important rôle in further constraining the sterile neutrino degrees of freedom. For completeness, we also display in Fig. 14 the LNV decays of D mesons into pseudoscalar and vector meson final states (accompanied by leptons of same or different flavour), D → ℓ α ℓ β V , with V = ρ, K * .

K meson decays
We conclude by presenting in Fig. 15 the prospects for the LNV decays of the K meson into a pion and a pair of leptons: for the ee and µµ final states the maximal expected branching ratios exceed the current experimental bounds over the majority of the considered sterile neutrino mass range (furthermore, the eµ channel is strongly constrained by the corresponding cLFV bound), although realistic observable rates are strongly suppressed by the upper bound on the sterile neutrino lifetime. Kaons are nevertheless an excellent framework to probe new physics via  contributions to LNV decays, since the already existing stringent bounds reported in Table 2 are likely to be improved in the near future by the NA62 collaboration, with an expected sensitivity of O(10 −11 ) for the ee and eµ channels, and O(10 −12 ) for the µµ channel 11 : should any hint of LNV be manifest in kaon decays, the results here presented can be useful in disentangling between the sterile neutrino hypothesis or other different mechanisms at its origin. Figure 15: Branching ratios of LNV decay modes of K mesons, K + → ℓ + α ℓ + β π − , as a function of the (mostly sterile) heavy neutrino mass, m 4 . Pale blue, yellow and green curves (surfaces) respectively denote the maximal (allowed) values of the BR(K + → ℓ + α ℓ + β π − ), with (ℓ + α ℓ + β ) = (e + e + , µ + µ + , e + µ + ); coloured dashed curves denote the corresponding maximal values of the BRs once cτ 4 ≤ 10 m is imposed. The coloured horizontal lines correspond to the present experimental bounds (cf. Table 2).

Tau-lepton decays
Due to its large mass, a tau lepton can also decay semileptonically; the kinematically allowed channels (always for the case of a primary production of an on-shell mostly sterile heavy neutrino) comprise several ∆L = 2 final states, including electrons and muons, as well as light kaons and pions.
The two panels of Fig. 16 display the branching ratios for final states with either an electron or a muon. It is interesting to notice that for the decay mode τ − → µ + π − π − , a heavy neutrino with a mass between 0.4 GeV and 0.55 GeV can be at the origin of decay widths already in conflict with experimental observation. Lepton number violation, in association with the τ − µ sector, thus emerges as a new constraint that must be taken into account in current analyses. Given the sizeable branching fractions for both µ + π − π − and µ + π − K − final states, any improvement in the experimental sensitivity (at LHCb, or Belle II) will render the LNV and cLFV observables a source of stringent constraints for models with additional Majorana neutrinos with masses below 1 GeV. Conversely, these channels might also play a discovery rôle for new physics sources of LNV.
3.2 Impact for the effective mass matrix: constraining m αβ ν One of the most important implications of the study of semileptonic LNV decays is the extraction of bounds on the relevant entries of the 3 × 3 effective neutrino Majorana mass matrix, m αβ ν . With the exception of the decays leading to ∆L e = 2, in which case the most constraining bound on m ee ν stems from neutrinoless double beta decays, the results we present emphasise to which extent LNV semileptonic decays allow to constrain the effective neutrino mass matrix.
As mentioned in the Introduction, past analyses have contributed to the derivation of upper bounds for certain entries (for instance in [10,28,30,[37][38][39][40]); it is important to notice that contrary to the neutrinoless double beta decay process, the dependence on the exchanged momentum (see Eq. (7), Section 2.4) -and hence on the masses of the different states (fermions and mesons) -does not favour a generalisation of the bounds. In this section, we thus discuss the expected sensitivity ; coloured dashed curves denote the corresponding maximal values of the BRs once cτ 4 ≤ 10 m is imposed. The coloured horizontal lines correspond to the present experimental bounds (cf. Table 5).
to the entries of the effective neutrino mass matrix, as inferred from the distinct three-body LNV semileptonic decays here addressed.
We begin by discussing the allowed values of the effective neutrino mass as a function of the mass of the heavy (mostly sterile) fermion as obtained from specific semileptonic decays. In particular, we choose as illustrative examples the predictions for the six independent entries of m αβ ν obtained from the LNV decays B → ℓ α ℓ β π (allowing for a large interval of sterile neutrino masses). These are complemented by the constraints for m µµ ν extracted from D → µµπ LNV decays, and m τ µ ν , which is directly inferred from semileptonic tau LNV decays (τ → µππ). Depending on the heavy neutrino length of flight (which is in turn a consequence of its mass and mixings -see Section 2.2.2), the LNV decays might occur outside a detector of finite size; leading to the results of all subsequent figures, we systematically exclude regimes which would be associated with L flight ν 4 10 m (denoted also by grey points). In Fig. 17, the different panels display the values of |m αβ ν | for the m 4 intervals kinematically allowed. (Leading to this figure, and to all subsequent ones, the underlying scan is similar to the one described for Fig. 12 in Section 3.1.) Two regimes can be identified for the experimentally allowed points (blue): the contributions arising from the light (active) neutrinos, associated with the saturation (mostly independent of m 4 ) observed for the smallest values of |m αβ ν |; the contributions from the heavy additional state, which span over several orders of magnitude (in GeV). The direct and indirect bounds on the branching ratios (as discussed in the previous subsection), typically translate into upper bounds of around 10 −3 to 10 −4 GeV for the entries of the effective mass matrix, with the exception of the τ τ entry. This already implies an improvement with respect to the results of [10,[36][37][38]. In particular, our bounds ameliorate those of [10,36] by 2 up to 7 orders of magnitude for the µµ, eτ and µτ entries 12 .
Other than being dependent on the sterile neutrino mass and mixings, notice that these bounds are intrinsically related to the actual LNV process. This is manifest when comparing the µµ panel of Fig. 17 with the left panel of Fig. 18, in which the bounds on |m µµ ν | are displayed -but now 12 Further bounds on LNV µ − e transitions, and hence on |m eµ ν | can be obtained from searches for lepton number violating neutrinoless µ − e conversion in Nuclei (µ − − e + , N); these were reported to be O(10 −2 GeV) by [36].
arising from decays of D mesons, D → µµπ: the bounds are considerably stronger, leading in general to |m µµ ν | 10 −6 GeV. Likewise, the bounds on |m µτ ν | obtained from LNV B decays can be compared with those obtained from semileptonic tau decays (right panel of Fig. 18): in this case both channels lead to similar constraints.
A final remark concerns the last entry, m τ τ ν , for which the available bounds are typically much less stringent. Although a same-sign tau pair in the final state is not always kinematically accessible, there are nevertheless decays, such as B → τ τ π, which allow to constrain m τ τ ν . As can be seen from the relevant panel of Fig. 17, one can already achieve important constraints, |m τ τ ν | 10 −2 GeV reflecting a significant improvement with respect to former results (bounds on |m τ τ ν | 10 4 GeV were suggested from an analysis of LNV conducted for searches at HERA [37]).
Complementary information (and an overview of the predictions for m αβ ν over the sterile neutrino parameter space) is displayed in Fig. 19: coloured isosurfaces for log(m αβ ν ) are identified in the parameter space generated by the heavy neutrino mass and the relevant combination of matrix elements, |U ℓα4 U ℓ β 4 |. Solid surfaces reflect excluded regimes due to the violation of experimental/observational bounds, already summarised in Fig. 3, or due to leading to L flight ν 4 10 m; the coloured lines denote the sensitivity reach of several future experiments (high energy and high intensity).
As can be seen from the different panels, the constraints on the mixings of the heavy neutrino to the active states can be sufficiently weak to allow sizeable contributions to the effective mass. In agreement with the discussion of Fig. 16, it is manifest that bounds arising from LNV tau decays already contribute to exclude small regions of the sterile neutrino parameter space (pink surfaces on bottom panels). The same analysis, but without excluding the contribution to LNV decays of long-lived heavy neutrinos is presented in Fig. 20. Since regimes associated with smaller mixings are now allowed, the surfaces corresponding to the regimes experimentally excluded due to conflict with meson and tau semileptonic bounds are also more important.
It is also interesting to notice that regimes associated with large effective masses (e.g. m µµ ν , as inferred from B → µµπ) are potentially within reach of future experiments, such as NA62, SHiP [115] and FCC-ee.

Discussion and concluding remarks
In this work we have revisited the prospects of minimal SM extensions via sterile fermions in what concerns the three-body LNV semileptonic decays of mesons and of tau leptons. Motivated by the recent progress in both experimental searches (new bounds on LNV decays as well as other relevant observables) and theoretical inputs (improved determination of decay constants, among others), we have evaluated the prospects of a large array of lepton number violating decay modes, leading to charged leptons of same or different flavours in the final state. Working in the framework of a toy-like extension of the SM via one sterile Majorana fermion, we have focused on the interesting (and promising) regime in which the latter state is on-shell. Our analysis further encompasses a full computation of the heavy fermion decay width, as well as the additional constraints of having the propagating neutrino decaying within a realistic (finite) detector.
After having carried an up-to-date evaluation of the sterile neutrino parameter space, including all available observational and experimental bounds, we have discussed to which extent the bounds on three-body LNV decays can further contribute to constrain the sterile neutrino mass and activesterile mixing angles. Our findings reveal that current data on semileptonic tau decays already allow to exclude certain regimes -in particular for heavy (mostly sterile) neutrino masses below 1 GeV. Lepton number violating K meson decay modes offer a remarkably rich laboratory: in addition to already excluding part of the sterile parameter space (based on existing bounds), these processes are expected to play a major constraining rôle with the future sensitivity of NA62. The results obtained for other decay modes (as for example D s → µµπ) further suggest that in the near future these decays could be within experimental sensitivity, and thus offer additional means to constrain the active-sterile mixing angles in the relevant sterile mass interval.
The study of the distinct LNV decay channels (leading to charged leptons of same or distinct flavour in the final state) has allowed to infer bounds on the different independent entries of the effective neutrino mass matrix. In our work we have used a general definition of m αβ ν (allowing to encompass the standard 0ν2β effective mass definition). The bounds obtained vary depending on the actual decaying particle and on the final states (strongly dependent on the mass and mixings of the Majorana mediator); although for the ee entry, neutrinoless double beta decay remains the most stringent bound, one can already derive bounds typically below (or around) O(10 −3 GeV). With the exception of the τ τ entry of the effective Majorana mass matrix -for which the derived bounds are weaker, m τ τ ν O(10 −2 GeV) -, the processes here considered constrain the remaining entries to lie below O(10 −3 GeV).
We conclude by presenting a global overview of the most relevant findings of our work, in particular the overall prospects of such a minimal SM extension regarding the LNV observables here addressed and the constraints on the Majorana effective neutrino mass matrix. In Fig. 21, we thus depict the semileptonic LNV branching ratios for different channels, as a function of the (corresponding) relevant entry of the effective neutrino mass matrix. If present, the horizontal lines denote the expected sensitivity of the corresponding experiment.
The near future experiments dedicated to neutrinoless double beta decays are expected to bring down the sensitivity to about |m ee ν | 0.01 eV [60]. Although the future sensitivity of LHCb, Belle II or NA62 might allow to improve upon the current constraints, lepton number violating semileptonic meson decays are unlikely to improve upon the 0ν2β bound for m ee ν . The two lower panels of Fig. 21 offer a critical comparison of the prospects of very distinct semileptonic decays regarding constraints on m αβ ν , more precisely on m µτ ν : despite the comparatively smaller kinematically allowed phase space for the τ decays 13 (which also accounts for the lower density 13 Notice that while the first three panels of Fig. 21 are associated with the same colour scheme for the m4 regimes, . The colour scheme of the in-laid points denotes regimes for the mass of the heavy (mostly sterile) neutrino, m 4 in GeV; grey points denote exclusions due to the violation of experimental/observational bounds, as well as having an excessively long-lived sterile state. Green crosses denote points which are excluded due to excessive contributions to LNV decays.
In this study we have considered a simple toy-like model, in which the SM is minimally extended by a massive Majorana sterile fermion, which is assumed to have non-negligible mixings with the active states (no hypothesis made on the mechanism of neutrino mass generation). The results here obtained -both regarding the updated constraints on the parameter space and the predictions regarding the LNV branching fractions and reconstructed effective mass -can nevertheless be interpreted as "benchmark" results for complete models of neutrino mass generation in which the SM is extended via n s sterile neutrinos (as low-scale type I seesaw models, and variants). Notice however that constructions in which the smallness of neutrino masses is explained via the smallness of a lepton number violating parameter and/or in which the heavy states form pseudo-Dirac pairs, lead to a suppression of the expected LNV rates. For instance, this is the case of the Inverse Seesaw [116], Linear Seesaw [117], specific realisations of the νMSM [118] and low-scale realisations of the type-I seesaw mechanism [119]. Further observables, and means to constrain other relevant parameters, might also be explored in association with the three-body semileptonic LNV decays here addressed, but this lies beyond the scope of the current work.
As shown in our study, such a minimal SM extension already leads to contributions to LNV observables (involving leptons of same or different flavour) which are close -if not in conflict!with current data. With the improvements of the experimental sensitivities, it is possible that one of these decays will be measured; in addition to confirming the Majorana nature of the mediator, such an observation would allow to test simple SM extensions via sterile neutrinos.

A Extensions of the SM via n S sterile fermions A.1 Formalism
Due to possible mixings with the light (mostly active) neutrinos, extensions of the SM via sterile states open the door to the violation of lepton flavour in both neutral and charged leptonic currents [120,121]. After electroweak symmetry breaking, and in the charged lepton's physical basis, the addition of n S sterile (Majorana) neutrinos leads to the following modification of vector and scalar currents: , In the above, g w denotes the weak coupling constant, cos 2 θ w = M 2 W /M 2 Z , P L,R = (1 ∓ γ 5 )/2, and m i are the physical neutrino masses (light and heavy); the indices α denote the flavour of the charged leptons, while i, j = 1, . . . , 3 + n S correspond to the physical (massive) neutrino states. The (SM) vector and axial-vector Z-couplings of charged leptons are parametrised by the C V and C A coefficients, defined as C V = −1 + 4 sin 2 θ w and C A = −1. The rectangular 3 × (3 + n S ) mixing matrix, U αj , encodes the mixing in charged current interactions (corresponding to the (unitary) PMNS matrix, U PMNS in the case of n S = 0); its upper 3 × 3 block, usually denoted U PMNS , denotes the mixing between the left-handed leptons. Finally, the matrix C parametrises lepton flavour violation in neutral currents, We notice that in addition to the vector and scalar currents above referred to, the interactions with neutral and charged Goldstone bosons are also modified: The modification of neutral and charged lepton currents as a consequence of the addition of sterile neutrinos to the SM content, opens the door to new contributions to a vast array of observables, possibly in conflict with current data. These are summarised in Appendix A.3, and will be applied throughout our phenomenological analysis.
A.2 Theoretical frameworks -the simple "3+1 model" A number of theoretical frameworks -ranging from minimal extensions to complete theoretical constructions -call upon sterile fermions: the latter are present in several mechanisms of neutrino mass generation, which in addition to accommodating neutrino data, also address in the baryon asymmetry of the Universe and/or put forward a viable dark matter candidate; likewise, sterile fermions can be minimally added to the SM in a simplified, toy-model approach to evaluate their impact regarding flavour and/or lepton number violation. In our analysis we will consider the latter approach, which we proceed to describe.
The SM field content is enlarged via the addition of a single Majorana sterile neutral fermion which, as noticed before, can also be interpreted as encoding the effects of a larger number of states possibly present in the underlying complete new physics model; this first phenomenological approach -minimal "toy model" -relies in an ad-hoc construction, which makes no assumption on the mechanism of neutrino mass generation, thus allowing to decouple the generation of neutrino masses (possibly occurring at higher scales, or even arise from interactions not calling upon the lighter sterile state) from the mechanism responsible for flavour violation and lepton number violation at low-energies.
As done in other works (see, e.g. [111][112][113][114]), in the present study we will rely on a simple toy model which is built under the single hypothesis that interaction and physical neutrino eigenstates are related via a 4×4 unitary mixing matrix, U ij . Other than the masses of the three light (mostly active) neutrinos, and their mixing parameters, the simple "3+1 model" is parametrised by the heavier (mostly sterile) neutrino mass m 4 , three active-sterile mixing angles as well as three new CP violating phases (two Dirac and one Majorana).

A.3 Constraints on sterile fermions
Due to the presence of the additional sterile states, the modified neutral and charged lepton currents might lead to new contributions to a vast array of observables, possibly in conflict with current data. These SM extensions via sterile fermions must be then confronted to all available constraints arising from high-intensity, high-energy and cosmological observations. Sterile states, with a mass above the electroweak scale, can have sizeable decay widths, a consequence of being sufficiently heavy to decay into a W ± boson and a charged lepton, or into a light (active) neutrino and either a Z or a Higgs boson. One thus imposes the perturbative unitarity condition [122][123][124][125][126][127], . Noticing that the leading contribution to Γ ν i is due to the charged current term, one obtains the following bounds [122][123][124][125][126][127]: where α w = g 2 w /4π, and C ii is given in Eq. (12). However, this constraint is not very relevant in the present analysis, as we focus on mass ranges for which the sterile neutrino is produced on-shell from meson or tau decays.
Observational constraints on the sterile masses and their mixings with the active states arise from an extensive number of sources. If kinematically accessible, sterile neutrinos can be produced in laboratory experiments via the interactions in (11): the negative observation of these processes constraints the mixings with the electron [96,99,100], muon [97,99,100] and tau [98,100] flavours. Moreover, and other than requiring compatibility between the left-handed lepton mixing matrix U PMNS and the corresponding best-fit intervals 14 defined from ν-oscillation data [1][2][3][4][5][6][7][8], we also impose, when relevant, unitarity bounds as arising from non-standard neutrino interactions with matter, on the deviation ofŨ PMNS from unitarity [104][105][106]. Further constraints on the activesterile mixings (and on the mass regime of new states) arise from electroweak precision observables; these include new contributions to the invisible Z-decay width (addressed in [128][129][130][131]), which must comply with LEP results on Γ(Z → νν) [41]; moreover, any contribution to cLFV Z decay modes should not exceed the present uncertainty on the total Z width [41], Γ(Z → ℓ ∓ α ℓ ± β ) < δΓ tot . In our study we also take into account current limits on invisible Higgs decays (relevant for m νs < M H ), following the approach derived in [132][133][134]. Likewise, negative results from laboratory searches for monochromatic lines in the spectrum of muons from π ± → µ ± ν decays are also taken into account [11,135]. The new states (through the modified currents) induce potentially large contributions to cLFV observables; we evaluate the latter [109,110,113,121,[136][137][138][139][140] imposing available limits on a wide variety of observables (some of them collected in Tables 3 and 4). In addition to the cLFV decays and transitions, which can prove instrumental to test and disentangle these extensions of the SM, important constraints arise from rare leptonic and semileptonic decays of pseudoscalar mesons decays (including lepton universality violating, cLFV and lepton number violating modes); we include constraints from numerous K, D, D s , B modes (see [141,142] for kaon decays, [143,144] for D and D s decay rates, and [145,146] for B-meson observations), stressing that in the framework of the SM extended by sterile neutrinos particularly severe constraints arise from the violation of lepton universality in leptonic meson decays (parametrised by the observables ∆r K , ∆r π and ∆r τ ) [131,147]. (Due to being associated with less robust experimental bounds, we do not include in our constraints the observables R eτ π and R µτ π .) Finally, we also take into account the recent constraints on neutrinoless double beta decay [59]: should the sterile states be Majorana fermions, they can potentially contribute to the effective mass m ee [148], which we evaluate following [107,149].
A number of cosmological observations [135,[150][151][152] put severe constraints on sterile neutrinos with a mass below the GeV (in particular below 200 MeV). In our study we will in general explore regimes associated with heavier sterile states (m νs 0.5 GeV) so that these constraints are not expected to play a relevant rôle.

B Computation of lepton number violating decay widths
In this Appendix we describe the theoretical computation of the LNV decay widths addressed in our work, in particular semileptonic tau decays and semileptonic decays of mesons into pseudoscalar and vector meson three-body final states.

B.1 Semileptonic tau-lepton LNV decay widths
Consider the LNV tau decay into two mesons and a charged lepton, mediated by the sterile neutrino ν 4 , and where ℓ = e or µ. The amplitude for this process can be computed as in which V M i and f M i are the CKM matrix and decay constant corresponding to the final state meson M i , U ℓ4 is the active-sterile mixing, m 2 ij are the momentum variables defined by m ij ≡ (k i + k j ) 2 and Γ 4 is the total decay width of the sterile neutrino ν 4 (the computation of the latter is described in Appendix C). Defining the first and second terms in Eq. (16) as iM τ 1 and iM τ 2 , the squared amplitude, spin-averaged over the initial state and spin-summed over the final state, is given by in which each term is given by Since we are interested in regimes close to a resonance (m 2 31 ≈ m 2 4 or m 2 23 ≈ m 2 4 ), the narrow width approximation can be applied as a good approximation. In this case, the propagator in the amplitude can be replaced by a δ-function as The interference term of Eq. (20) contains in fact two resonances; these can be split into two separate parts (f 1 and f 2 ), each including only one resonance 15 , as with f τ i defined by This allows to remove one of the integrals in Eq. (21) by the δ-function introduced in Eq. (23). After applying the narrow width approximation, the remaining integration intervals are given by with E * 1 = (m 2 τ − m 2 4 − m 2 M 1 )/(2m 4 ), and E * 3 = (m 2 4 − m 2 M 2 + m 2 ℓ )/(2m 4 ). Notice that the mass range of the sterile neutrino ν 4 is limited to respectively for the integrals f τ 1 and f τ 2 , as a consequence of the narrow width approximation.

B.2 Widths of semileptonic LNV meson decays
Let us now consider the decay where both M 1 and M 2 are pseudoscalar mesons. The corresponding amplitude is computed as 15 Single-Diagram-Enhanced multi-channel integration [153].

D Computation of widths for charged meson cLFV semileptonic decays
While in SM extensions via additional sterile fermions the cLFV semileptonic decay of neutral mesons, i.e. M 0 1 → ℓ ± α ℓ ∓ β M 0 2 , necessarily occurs at loop-level, the analogous decays of charged mesons can proceed via the so-called "double W " diagrams, already at the tree-level.
Should their mass allow for the production and propagation of real, on-shell neutrinos on the s-channel (as is the case for the regimes we are interested in our study) then one can write which is valid in the limit of narrow-width approximation for the neutral fermion; the computation of the different heavy neutrino decay widths leading to Γ tot ν k has been detailed in Appendix C. The width of a leptonic charged meson decay can be cast as (see, for example [131]) in which the kinematical function λ(1, x ℓ , x ν k ) is given in Eq. (44), and all other quantities have been previously defined. For completeness, we also include here the full charged meson decay width, for the generic propagation of virtual, massive neutrinos [15]: where F(s, x M 2 , x ℓα , x ℓ β , x ν k ) = 1 s λ 1/2 (s, x M 2 , x ℓα ) λ 1/2 (s, 1,