Sensitivity of Lepton Number Violating Meson Decays in Different Experiments

We study the discovery prospect of different three body lepton number violating~(LNV) meson decays $M_{1}^{-}\to\ell_{1}^{-}\ell_{2}^{-}M_{2}^{+}$ in the framework of right handed~(RH) neutrino extended Standard Model~(SM). We consider a number of ongoing experiments, such as, NA62 and LHCb at CERN, Belle II at SuperKEK, as well as at the proposed future experiments, SHiP, MATHUSLA and FCC-ee. The RH Majorana neutrino $N$ mediating these meson decays provides a resonant enhancement of the rates, if the mass of $N$ lies in the range $(100\, \text{MeV}-6\, \text{GeV})$. We consider the effect of parent mesons velocity, as well as, the effect of finite detector size. Using the expected upper limits on the number of events for the LNV decay modes, $M_{1}^{-} \to\ell_1^{-}\ell_2^{-}\pi^{+}$~($M_{1}=B, B_c,D, D_{s}\,\text{and}\,K$), we analyze the sensitivity reach of the mixing angles $|V_{e N}|^{2}$, $|V_{\mu N}|^{2}$, $|V_{\tau N}|^{2}$, $|V_{e N}V_{\mu N}|$, $|V_{e N}V_{\tau N}|$ and $|V_{\mu N}V_{\tau N}|$ as a function of heavy neutrino mass $M_{N}$. We show that, inclusion of parent meson velocity can account to a large difference for active-sterile mixing, specially for $D$, $D_s$ meson decay at SHiP and $K$ meson decay at NA62. Taking into account the velocity of the $D_s$ meson, the future beam dump experiment SHiP can probe $|V_{eN}|^2 \sim 10^{-9}$. For RH neutrino mass in between 2 - 5 GeV, MATHUSLA can provide best sensitivity reach of active-sterile mixings.


INTRODUCTION
The discovery of neutrino oscillations [1] in a series of oscillation experiments have confirmed that neutrinos have non-zero masses and non-zero mixings. The solar and atmospheric mass splittings are O(10 −5 ) and O(10 −3 ) eV 2 , while the three mixings are θ 12 ∼ 33 • , θ 23 ∼ 45 • and θ 13 ∼ 9 • . These observations indicate, at least two of the three SM neutrinos have non-zero masses. The absolute scale of the neutrino masses are yet unknown.
The sum of masses of three active neutrinos are bounded from cosmological observation as i m ν i < 0.23 eV [2]. One of the most natural ansatz to explain small neutrino masses is the seesaw mechanism, where the dimension-5 operator [3] with lepton and Higgs doublets generates Majorana mass of light neutrinos through electroweak symmetry breaking.
This operator breaks global lepton number symmetry of the SM by two units. The other possibility is to generate Dirac mass terms for the SM neutrinos by including gauge singlet RH neutrinos in the theory. However, to explain eV mass of the neutrinos, this requires unnatural fine-tuning of Yukawa coupling to a very small value Y ν ∼ O(10 −11 ). The seesaw mechanism on the other hand is most appealing, as the tiny Majorana mass of the light neutrinos are inversely proportional to the cut-off scale of the theory. This large cut-off scale naturally explains eV mass of neutrinos. Seesaw can be realised in different beyond standard model (BSM) extensions, such as Type-I [4][5][6][7]/Inverse seesaw [8,9] with gauge singlet Majorana neutrinos, Type-II [10][11][12][13] seesaw with Higgs triplets, and Type-III seesaw [14][15][16][17] with fermionic triplet. For Type-I/Inverse seesaw, RH neutrinos can have Majorana/Quasi-Dirac masses, that can vary in wide ranges, starting from GUT scale down to GeV. The low scale seesaw models, that inherits lighter RH neutrino states have higher discovery prospect, as they can be tested in a wide range of experiments.
Another interesting probe for Majorana neutrinos of hundreds of MeV-few GeV masses and their mixings are the lepton number violating (LNV) rare meson decay processes, such as M − 1 → − of the ongoing and future experiments. In secs. 9, we present our combined limit from all the considered meson decays and give the comparison with other existing constraints on the mixing angles. Finally in sec. 10, we provide our conclusions. In the Appendix, we provide details of the RH neutrino decay width calculations.

THE MODEL
We extend the SM to include additional RH neutrinos N . The heavy neutrinos can generate light neutrino masses through seesaw. For simplicity, we consider only one RH neutrino and carry out our analysis. The mixing of N with the active neutrinos are given by the following expression, where ν m and N m are the mass eigenstates. We denote the mixing between the standard flavour neutrino ν ( = e, µ, τ ) and the heavy mass eigenstate N by V N . Due to the mixing, the charged and neutral currents in the lepton sector gets modified and can be written as For the purpose of phenomenology, we consider the mass and mixings of N as free parameters, constrained only by experimental observations. Note that, adding only one RH neutrino is not enough to correctly reproduce the neutrino oscillations parameters, namely two mass square differences and the mixings. In our considered model, we can add two more RH neutrinos to generate the neutrino masses and consider two of them to be heavy enough such that only one RH neutrino lies in the mass range 0.140 GeV ≤ M N ≤ 6 GeV.

PROCESS
The RH neutrino N , if a Majorana state can mediate the LNV process M − 1 → − 1 − 2 M + 2 . The Feynman diagrams for these decays are shown in Figs. 1 and 2. The diagram in Fig. 2 will give a very small contribution, as this is not a resonance production diagram. Note The Feynman diagrams for the lepton number violating meson decays. These processes can produce resonance enhancement. See text for details that the diagrams with light neutrino exchange are also present but the contributions will be negligibly small as they will not be resonantly enhanced. The decay amplitude for the Fig.1 can be expressed as 188 GeV and f Bc = 0.436 GeV [64]. M N , Γ N are the mass and decay width of the heavy neutrino N . Finally, the total decay rate is given by In Eq. 5, n = 2 for identical final leptons, otherwise n = 1 and d 3 (PS) is the three body phase space. the following channels contribute to the total decay width of heavy neutrino: • N → − P + , where = e, µ, τ , and P + = π + , K + , D + , D + s , B + (for = e, µ).
Hence, the total decay width is given by As the RH neutrino is Majorana, the charge conjugate processes are also allowed and the decay rate is same, hence the 2 factor is included for some of the channels. We can parameterize the above decay width as follows where, a e , a µ and a τ are functions of the Majorana neutrino mass and hence will differ from mode to mode. We show the total decay width of heavy neutrino N in the left panel of Fig. 3, for the choice of mixing |V eN | 2 = |V µN | 2 = |V τ N | 2 = 1. Even with mixing angle equal to 1, the decay width lies in the range 10 −16 GeV ≤ Γ N ≤ 10 −7 GeV for the M N mass range, 0.140 GeV ≤ M N ≤ 6 GeV. Hence, we can safely use the narrow-width approximation, . We also show in the right panel of Fig. 3, the different coefficients

PARENT MESON VELOCITY AND FINITE DETECTOR SIZE EFFECT
For the mass range 0.140 GeV ≤ M N ≤ 6 GeV, the RH neutrino produced in these LNV meson decays are on shell. The RH neutrino produced in meson decays M 1 → 1 N , propagates and decays after traveling some distance from the production point. This is the decay length L N of the RH neutrino N and it depends on the total decay width of N . If L N is greater than the actual size of the detector, then N will decay outside the detector and the signature M 1 → 1 2 π cannot be observed. For a particular experiment the detector size is finite. Hence, when calculating the signal events, we need to take into account this finite size detector effect by the probability factor P N , of N to decay within the detector.
In general, this probability factor can be written as where L N = p N M N Γ N , L D is the detector length, p N is the three momentum of N . Defining x = L D L N , it is obvious that for a very large detector length L D and small decay length L N , P N = 1 − exp(−x) ≈ 1. Note that the probability factor depends on three momentum p N , which in turn depends on the velocity of decaying meson M 1 . Hence to incorporate the probability factor correctly, we need to use the correct velocity of the parent meson M 1 in each of the experiments. If the parent meson M 1 decays at rest, three momentum p N is fixed and is given by . For the case of parent meson M 1 produced with fixed boost β, the energy of N is then given by, where E * N is the energy of N in rest frame of M 1 which is given as and θ * N is the emission angle of particle N in the rest frame of M 1 , which is measured from the boost direction β. The energy E N of the N in the boosted M 1 frame lies within the range, Similarly we can derive the range of p N ∈ [p − N , p + N ] from Eq. 9 using the relation p ± In this section, we show how x, L N depends on parent meson velocity and compare to the case of parent meson decay at rest. For the case of meson decay at rest, p N = p * N and for meson decay with non-zero momentum p M 1 , we take p N = compare. With the assumption of |V eN | 2 =|V µN | 2 =|V τ N | 2 , we can write the decay length L N and x as In Figs 4 and 5, we have shown the variations of L N .|V N | 2 and x |V N | 2 as a function of RH neutrino mass M N in B, D s and K meson decays. To do the comparison with the meson decay at rest (p M 1 = 0), we take p B = 45 GeV (FCC-ee) [65], 58 GeV (SHiP) [66]; p Ds = 58 GeV (SHiP) and p K = 75 GeV (NA62) [67]. From these two figures it is clear that decay length increases (hence x decreases) for fixed mixing angle in the case of meson decays in flight compared to meson decay at rest. Hence, the probability of RH neutrino P N to decay inside the detector is smaller in the case of meson decay in flight comapare to meson decay at rest. As a result, compare to meson decay at rest, in the case of meson decay in flight we get a rather loose bound on mixing angle from the expected signal events.

SIGNAL EVENTS
The sensitivity reach for the LNV decay modes in a particular experiment depends on the number of the parent mesons M 1 's produced (N M − 1 ), their momentum ( p M 1 ) and the branching ratio for these mesons to the LNV modes. Assuming the parent meson M 1 decays at rest, the expected number of signal events is [68]: the factor 2 is due to inclusion of the charge conjugate process M + 1 → + 1 N and P N is the detector probability which is given by For the case of meson decay in flight the RH neutrino energy E N lies in range according to Eq. 9 and follows a flat distribution as: Hence to calculate the total number of events for M − 1 → − 1 − 2 M + 2 in the lab-frame we need to integrate within the range of E N as where is the detector probability after taking into account the parent meson M 1 velocity.
Since the LNV meson decay rates will be extremely small, the expected number of signal events for these processes can be assumed to follow a Poisson distribution. Following Ref. [69] and assuming zero background events, we derive the average upper limit on the number of events at 95% C.L., assuming number of signal events to be N event = 3.09.
Note that the number of events given in Eqs. 12 or 13 are functions of the mass parameters M N and mixing V N . Equating the numerical upper limit on the number of events to the theoretical expressions, we get constraints on mixing angle V N , corresponding to specific when deriving these bounds using Eq. 12 and 13. LHCb are smaller than the number of B mesons, this mode being less suppressed with respect to B + → + 1 + 2 π − , gives tighter constraints on the mixing angles. The produced mesons will decay in flight, carrying a momentum of order of 100 GeV in forward direction [76]. We take the detector length L D ≈ 20 m.

NA62
NA62 is an ongoing experiment at CERN that will produce a large number of K + mesons [67]. The primary SPS 400 GeV proton beam, aims on a target, produce a secondary high intensity hadron beam with an optimum content of K + (≈ 6%). The expected number of K + decays in the fiducial volume is 4.5 × 10 12 per year. Assuming three years of running, N K + = 1.35 × 10 13 . The detector length L D ≈ 65 m and the produced K + mesons will decay in flight, carrying a momentum of 75 GeV.

SHiP
The SHiP experiment is a newly proposed general purpose fixed target facility at the CERN SPS accelerator [77]. A 400 GeV proton beam will be dumped on a heavy target for the duration of five years. One of the primary goal of the experiment is to use decays of charmed mesons to search for heavy sterile neutrinos using the decay mode D + s /D + → + + π − . One can easily estimate the number of charmed meson pairs that are expected to be produced in this experiment as [78], where X cc is the cc production rate, N P OT = 2 × 10 20 is the number of proton-target interaction. The relative abundances R of charmed mesons, such as, D and D s are 30% and 8%, respectively. Hence, the expected number of D and D s mesons are N D + = 1.02 × 10 17 and N D + s = 2.72 × 10 16 , respectively. This very high intensity of the charmed mesons will permit to set tight constraints on mixing angle at SHiP. There will also be large number of RH neutrinos which are coming from the meson decays or W and Z boson decays has very large decay length and this has been already studied in [80]. For the meson decay case, they have considered the decay modes B → D N , B → N and D → K N and after including the probability of RH neutrinos to decay visibly within the MATHUSLA detector, they derived the constraints on mixing angles. In this study, we have considered the meson decays B → 1 2 π and D → 1 2 π for MATHUSLA. For the number of B and D meson productions we followed the Ref. [81]. The result of their detailed simulation suggests that number of B and D meson production within the geometric acceptance of the MATHUSLA detector are 5.7 × 10 14 and 5.4 × 10 13 , respectively. Their simulation also gives the average γ factor of the B and D mesons as γ B = 2.3 and γ D = 2.6 from which we can derive the average momentum of the mesons. The detector length is taken to be 38m.

Belle II
The asymmetric SuperKEKB facility is designed to collide electron and positron beams such that the centre of mass energy is in the region of the Υ resonances. An upgrade of Belle, the newly completed Belle II detector is expected to collect data samples corresponding to an integrated luminosity of 50 ab −1 by the end of 2024 [82]. The expected number of charged mesons will also be accessible, with N D + = 3.4 × 10 10 and N D + s = 10 10 [84]. A direct search for heavy Majorana neutrinos in B-meson decays was performed by Belle collaboration using a data sample that contained 772 × 10 6 BB pairs (at 711 fb −1 ) [85]. At KEKB as well as superKEKB, the energies of the e + , e − beams are sufficiently low so that the momentum of the produced B mesons as well as that for the charmed mesons will not be appreciable and the suppression from high momentum of the decaying mesons in the number of events will be absent.
where N Z ∼ 10 13 , Br Z → bb = 0.1512 [86], f u = 0.410 [87] is the fraction of B + from b quark in Z decay. The B mesons produced at FCC-ee will have an energy distribution peaked at E B + = M Z 2 . Hence we can calculate the number of signal events using Eq. (13), where the detector length is taken to be L D = 2 m.

RESULTS
In Fig. 6, we show how the velocity of the parent mesons affect the sensitivity reach of the mixing angles. We consider a number of ongoing and future experiments, such as, FCC-ee, SHiP to explore B → eeπ, SHiP for D s → eeπ, NA62 for K → eeπ, and LHCb, SHiP for B c → τ τ π meson decays. To derive the bounds/sensitivity on the mixing angle as a function of RH neutrino mass M N , we use Eq. 12, and Eq. 13, for meson decay at rest and in flight, respectively. For all of the above decays, the obtained bounds/future sensitivity on the mixing angles are rather loose in case of meson decays in flight compared to meson decays at rest. As an example, for the case of B (FCC-ee, SHiP) and D s (SHiP) meson decays, there is approximately one order of magnitude difference between the two results. The result for K meson decay at NA62 differs by two order of magnitude. Hence, of charmed meson productions, the future experiment SHiP will be able to probe |V eN | 2 , |V µN | 2 ∼ O(10 −9 ) and |V eN V µN | ∼ O(10 −9 ) in the mass range 0.14 GeV < M N < 1.9 GeV. at LHCb, respectively. Note that, B, B c → τ τ π meson decays constraint the mixing angle |V τ N | 2 in the mass range, where it has so far been unconstrained by any of the τ or other meson decays. In spite of larger number of D production, compared to D s meson at SHiP, the suppression from the weak annihilation vertex in the case of D s meson V cs is less compared to D meson V cd . As a result of this, tightest limits on the mixing angles will be provided by the D s meson decays in the relatively lower mass range.
Note that, if both the like sign di-leptons are not of the same flavor ( 1 = 2 ), then the process is not only lepton number violating, but also lepton flavor violating. Further, if the distance between N production and decay points is large enough, then the two processes, M 1 → 1 N followed by N → 2 π will be separated. Assuming this separation, the two processes M 1 → 1 2 π and M 1 → 2 1 π can be distinguished. While deriving the bounds on the mixing angles for the case of 1 = 2 , we are assuming this separation in our study.
This is justified as the decay width Γ N is very small (hence the lifetime is very large) in the mass range of interest. The allowed mass range of N for the decay modes M 1 → 1 2 π respectively. We consider both of the channels M 1 → 1 2 π and M 1 → 2 1 π to derive the bound on the mixing angle |V 1 N V 2 N |.
One important point to note is that we have considered idealized detector with 100% detection, reconstruction efficiencies etc to derive the constraints on the mixing angles. The realistic constraints are expected to be weaker and will only be feasible through searches by the experimental collaborations, incorporating the detection, reconstruction efficiencies in actual experiment.

COMBINED SENSITIVITY REACH FROM MESON DECAYS AND COM-PARISON WITH EXISTING CONSTRAINTS
In this section, we discuss the future sensitivity reach from LNV three body meson decays. The combined limits represent the strongest limits obtained in different mass ranges of N . In Fig. 11, we show the combined sensitivity reach for |V eN | 2 by dark blue solid line. This corresponds to the tightest constraints that can be obtained from D s → eeπ mode (by SHiP) in the lower mass range 0.14 GeV < M N < 2 GeV, and from B → eeπ (by  Note that, for the lower mass range, SHiP can probe |V eN | 2 ∼ 10 −9 , while for higher mass range, MATHUSLA and LHCb can probe |V eN | 2 ∼ 10 −7 . In particular, the very near future accumulation of data (L = 300 fb −1 ) in LHCb can probe |V eN | 2 ∼ 10 −7 , around M N ∼ 5 GeV. The sensitivity reach of |V µN | 2 , as shown in Fig. 12 is very similar. The combined limit represents the constraint from D s → µµπ, B → µµπ and B c → µµπ decay modes, that can again be probed in SHiP, MATHUSLA and LHCb. For |V τ N | 2 , the best sensitivity reach |V τ N | 2 ∼ 10 −7 can be provided by MATHUSLA in B → τ τ π mode, while SHiP and LHCb can give similar sensitivity reach with the mode B c → τ τ π. The combined sensitivity reach has been shown in Fig. 13.
The future sensitivity of |V eN V µN |, |V µN V τ N | and |V eN V τ N | are shown in Fig. 14, Fig. 15 and Fig. 16, respectively. For |V eN V µN | mode, lower mass range up to M N ∼2 GeV can be probed by the channel D s → eµπ ( by SHiP) where sensitivity down to |V eN V µN | ∼ 10 −9 can be obtained. RH neutrino of higher mass M N ∼ 5 GeV and M N ∼ 6 GeV can be probed by B → eµπ mode (by MATHUSLA) and B c → eµπ mode (by LHCb), with sensitivity reach |V eN V µN | ∼ 10 −7 . For |V eN V τ N | and |V µN V τ N | mixings, the sensitivity for the active-sterile mixing angles are similar, as depicted in the Fig. 15 and Fig. 16. These can be probed in LHCb, MATHUSLA and SHiP. We note that, the future limits from LNV meson decays will be most sensitive in between 0.5 GeV < M N < 2 GeV for |V eN | 2 , |V µN | 2 , |V eN V µN |. For other mixing angles that involves τ in the final state, best limit can be obtained in relatively higher mass range M N ∼ 5 GeV. We stress that, LHCb, and future experiments -SHiP and MATHUSLA can probe mixing angle of the τ sector in a region, that is very loosely constrained.
The PS191 [115] limits for the electron and muon flavors are shown in Figs. 11 and 12, respectively. The JNIR [116] limit is represented by the black dashed line for the electron flavor in Fig. 11. The regions excluded by the present constraints are shaded in gray. The limits from GERDA [117] on the mass mixing plane in search of Majorana neutrinos from the neutrinoless double beta decay is represented by the dark cyan line in Fig. 11. Majorana heavy neutrino searches from the Meson decay in E949 [118] and NuTeV [119] can produce strong bounds on the heavy neutrino mass-mixing plane. Lepton-jet theoretical searches [120] for the Majorana neutrinos with muon flavor can also put strong bounds. These bounds are shown in Fig. 12 for the muon flavors. For the tau lepton, the bound in the corresponding mass region from EWPD [121][122][123] has been shown in Fig. 13. The decay of tau lepton into heavy Majorana neutrino and meson can also put bounds on the massmixing plane and can have prospective limits marked as B-factory [89,92] 2 . The NA62 [125][126][127] projection lines from the electron, muon and tau are shown in Figs. 11, 12 and 13, respectively. Such experiment can be performed in the kaon mode and beam dump mode [128]. This search is sensitive to the heavy neutrinos that are produced in weak decays [129,130] of mesons or tau leptons [66]. The upper limit on the mass mixing plane from We briefly summarize the current strongest experimental bounds on the mixing angles such as |V eN V µN |, |V eN V τ N |, and |V µN V τ N | for the Majorana heavy neutrinos in Fig. 14, 15 and 16, respectively for M N < 10 GeV. Strongest bounds from the BBN [66,106,107], PS191 [115], NuTeV [119], CHARM [108][109][110], DELPHI [111] are obtained from the Majorana heavy neutrino search for M N ≤ 10 GeV. The shaded region is excluded by the results obtained from these experiments. The prospective bounds from the µ → e (Ti) and µ → e (Al) are shown in Fig. 14 from [138]. The limits on the mixings from the τ decay into hadrons in association with electron and muon are shown in Figs. 15 and 16, respectively from BABAR [139]. The limits from the τ → eππ and τ → e − π + π + are shown in Fig. 15 and those obtained from τ → µππ and τ → µ − π + π + are shown in Fig. 16, respectively [140].
We stress that, in the relatively lower mass range, among the experimental constraints, the tightest constraint on the mixing angles |V eN | 2 and |V µN | 2 can be obtained from LNV meson decays. These are however still one order of magnitude weaker than the theory constraints from BBN and Seesaw. For relatively higher mass range, our combined bounds on the mixing angles |V eN V µN |, |V eN V τ N | and |V µN V τ N | are the tightest bounds. As we have discussed before, the LNV meson decays can probe the product of the mixings |V µN V τ N |, |V eN V τ N | in higher mass ranges M N ∼ 5 GeV, that are so far unconstrained. SHiP. Due to non-zero velocity of the parent meson, the probability of the generated RH neutrino to decay inside the detector changes. This alters the sensitivity reach by more than an order of magnitude for the above mentioned experiments.
We explore a number of channels, B/D/D s → µµπ, B/D/D s → eeπ, B c /D s → eµπ, B → τ τ π, B/B c → eτ π, B/B c → µτ π, and few others. We find that, for the mass range M N ∼ 1 GeV, future experiment SHiP can probe |V eN | 2 ∼ 10 −9 , while for mass range M N ∼ 5 GeV, future experiment MATHUSLA, and LHCb with 300 fb −1 integrated luminosity can probe |V eN | 2 ∼ 10 −7 . The sensitivity reach of |V µN | 2 is very similar to |V eN | 2 . For

APPENDIX
The different partial decay widths of the RH neutrinos N i are, In the above decay mode 1 = 2 .
where x i = m i M N with m i = m , m P 0 , m V 0 , m P + , m + V . The kinematical function are given by, Neutral current couplings of leptons are given by, Neutral current coupling of pseudoscalar and vector mesons are given by,