Sensitivity bounds on heavy neutrino mixing $|U_{\mu N}|^2$ and $|U_{\tau N}|^2$ from LHCb upgrade

Decays of heavy pseudoscalar mesons $B$, $B_c$, $B_s$ and $D_s$ at LHCb upgrade are considered, which produce either two equal sign muons or taus. In addition, we consider the analogous decays with opposite sign muons or taus. All these decays are considered to be mediated by a heavy on-shell neutrino $N$. Such decays of $B$ mesons, if not detected, will give in general stringent upper bounds on the heavy-light mixing parameter $|U_{\mu N}|^2$ as a function of the neutrino mass $M_N \sim 1$ GeV, principally due to the large expected number of produced mesons $B$. While some of the decays of the other mentioned mesons are attractive due to a weaker CKM-suppression, the expected produced number of such mesons is significantly smaller that that of $B$'s; therefore, the sensitivity bounds from such decays are in general comparable or less restrictive. When $\tau$ pairs are produced, only two types of such decays are significant: $B^{\pm}, B_{c}^{\pm} \to \tau^{\pm} \tau^{\pm} \pi^{\mp}$ (and $\tau^{\pm} \tau^{\mp} \pi^{\pm}$), giving us stringent upper bounds on $|U_{\tau N}|^2$; the other decays with a pair of $\tau$, such as $B^0 \to D^{(*)-} \tau^+ \tau^+ \pi^-$ (and $D^{(*)-} \tau^+ \tau^- \pi^+$), are prohibited or very suppressed by kinematics.


I. INTRODUCTION
Various theories which can explain the masses of the three light neutrinos suggest the existence of heavy neutrinos, often referred to as sterile neutrinos. The neutrinos in such theories are of the Majorana type, i.e., they are their own antiparticles. As a consequence, they can produce, apart from lepton number conserving (LNC) processes, lepton number violating (LNV) processes through on-shell mediation. Dirac neutrinos can participate only in LNC processes. However, at present there are many open questions in the neutrino sector, among them: (a) What is the nature of the neutrinos, i.e., are they Majorana or Dirac particles? (b) How many heavy neutrinos exist and what are the values of their masses? (c) What are the values of their heavy-light mixing parameters U N (where = e, µ, τ )?
In some of these processes, the neutrino mass can also be determined or constrained. If a considered pseudoscalar meson M decays into invisible channel, M → invisible, then it is difficult or impossible to extract or constrain the mass M N of the involved neutrino(s) N (via invisible width measurement), principally because we have competitive irreducible SM background M → νννν [17] 1 However, if a rare decay process of M is mediated by an on-shell neutrino N (M N ∼ 1 GeV) and, simultaneously, all or most of the final state particles are detectable, then the mass M N can be either determined (M 2 N = P 2 f ) or at least reasonably constrained. In this work we consider the rare decays M → N → π ± where = µ or τ , i.e., all the final state particles are charged and in principle detectable.
Furthermore, the phenomenon of oscillations [18] between the three light neutrinos has been observed [19]. If heavy neutrinos have almost degenerate mass, they can also oscillate [20]. A somewhat related phenomenon is the resonant CP violation which can arise in scattering processes [21], rare meson decays [8,[22][23][24] and rare τ decays [14]. There are several models with almost degenerate heavy neutrinos, among them the low-scale seesaw models [25] and the neutrino minimal standard model (νMSM) [26,27].
Various models with seesaw mechanism [28] and related models explain very low masses of the light neutrino sector, and contain heavy neutrinos with masses M N 1 TeV [29], M N ∼ 1 TeV [30] and M N ∼ 1 GeV [3,26,31]. In the latter type of models, not only the masses are relatively low, but the heavy-light mixing parameters |U N | 2 are often less suppressed than in the earlier seesaw models [29]. Our analysis will be made with a view to such scenarios, i.e., M N ∼ 1 GeV and |U N | 2 ∼ 10 −6 -10 −4 (for = µ, τ ).
In this work we continue and extend the analysis of Ref. [32] of the sensitivity bounds on the heavy-light mixing parameter |U µN | 2 in the LHCb upgrade, where the rare lepton number violating (LNV) decays of B-mesons with two equal sign muons (LNV) were considered. Here we recalculate these rare B-decays, with the updated input parameters [33], and calculate also the analogous decays with opposite sign muons (LNC) when the neutrino N is Dirac. If the intermediate neutrino N is Dirac, only LNC processes take place; on the other hand, if N is Majorana, both LNC and LNV processes take place. In addition, we consider in this work the rare LNV and LNC decays of the pseudoscalar mesons M = B c , B s and D s , with two muons; and the rare decays of M = B and B c with two taus. We consider that the decays are mediated by a heavy on-shell Majorana or Dirac neutrino N with mass M N ∼ 1 GeV. The considered decays with muons involve the meson decay into vector V or pseudoscalar meson S and an off-shell W * , where the latter gives a muon and and (on-shell) neutrino N : M → V (S)µN . The neutrino N propagates from the primary vertex and decays at the secondary vertex within the detector, producing a muon and a pion: M → V (S)µN → V (S)µµπ. The considered decays are: (a) for B mesons B → D * (D)µµπ, and B → µµπ; (d) and for D s mesons D s → Φ(η)µµπ, and D s → µµπ. When the two charged leptons are τ 's in the decays of the mentioned types, the kinematics allows only the consideration of two decays: B → τ τ π and B c → τ τ π.

II. THE FORMALISM USED
The analysis in Ref. [32], which we continue using here, was based on the formalism and expressions obtained in Ref. [34]. As in those references, we take the simplest representative scenario of heavy neutrinos, namely one heavy Majorana or Dirac neutrino N (with mass M N ∼ 1-10 GeV), which has suppressed (heavy-light) mixing parameters U N with the light flavor neutrinos ν ( = e, µ, τ ) Here, ν k (k = 1, 2, 3) are the three light neutrino mass eigenstates.
Here we summarize the expressions obtained in Ref. [34]. The decay widths for Here j denote charged leptons; later we will have 1 = 2 = µ or 1 = 2 = τ . The first factor, Γ (M → V (S) 1 N ), was calculated in Ref. [34] for the case of vector (V) and scalar (S) meson, and we refer to that reference for details. We note that the mass M N ∼ 1 GeV plays an important role. The second factor, Γ(N → 2 π), is well known, and is also given explicitly in Refs. [32,34]. The total decay width Γ N ≡ Γ(N → all) was given and evaluated numerically as a function of M N in Ref. [23] for the case of Majorana neutrinos, and in Ref. [8] for the cases of Dirac and Majorana neutrinos.
As pointed out in [32], the branching ratio of the process has to be multiplied by the probability P N of the intermediate N neutrino to decay within the detector. This factor has to be considered with care, because it depends on the kinematic parameter β N which is the velocity of N in the lab frame (Σ ). This β N depends in a nontrivial way on the quantities q 2 ,q , andp 1 , where q is the four-momentum of the virtual W * (→ 1 N ),q is its unitary direction vector in the M -rest frame (Σ ), andp 1 is the direction of of the produced 1 in the W * -rest frame (i.e., 1 N -pair rest frame, Σ) Here, E N is the energy of N neutrino in the laboratory frame (Σ ), and L is the effective length of the detector. 2 We refer to [32] for details. The true (effective) branching ratio Br eff of the process is then evaluated by integrating the corresponding differential decay rates multiplied by the decay probability P N and divided by Γ M The values of the total decay widths Γ M were taken from [33]: Γ B + = 4.018 × 10 −13 GeV; Γ B 0 = 4.330 × 10 −13 GeV; Γ Bc = 1.298 × 10 −12 GeV; Γ Bs = 4.326 × 10 −13 GeV; Γ Ds = 1.306 × 10 −12 GeV. The decay widths Γ N ≡ Γ(N → all) of N Majorana and Dirac neutrinos are obtained using the expressions of Ref. [23] (Appendix B and Fig. 2 there); cf. also [8] (Appendix A.3 and Fig. 2 there). They have the form where the factor K in Eq. (5) has all the dependence on the heavy neutrino mixing factors The (dimensionless) numbers N N are functions of the mass M N , N N (M N ) ∼ 1-10. They are determined by the formulas given in Appendix B of Ref. [23] which are based on the approach of Ref. [6]. The results for N N (M N ) are presented, e.g., as a curve in Fig. 2 of Ref. [23]. When the produced meson is vector (V ), the differential decay width dΓ(M → V 1 N )/(dq 2 dΩq dΩp 1 ) is complicated [34], due to both the vector nature of the produced meson V and the nonzero mass of N (and of 1 ); but it does not depend on the electric charge of 1 . The available literature does not shed light unequivocally on the latter point, for which we refer to the Appendix.
If no mesons V or S are produced, the differential decay width for M → 1 N is simpler, it depends only on the directionp N of the on-shell N in the M -rest frame; and since M is (pseudo)scalar, we have dΓ( . Further, the decay probability P N also depends only onp N . This then gives (for details, cf. ( [32,34]) III. NUMERICAL RESULTS FOR THE ACHIEVABLE UPPER BOUND LIMITS OF |U N | 2 ( = µ, τ ) AT LHCB UPGRADE We will consider the mentioned rare decays of M = B, B s , B c and D s at LHCb upgrade. We take into account that the produced B mesons have a specific distribution of the momentum (| p B |) lab ≡ | p B | in laboratory (Σ ) frame, as given in Fig. 1(a). 3 In practice, we separated this distribution in ten bins with equal weight (i.e., equal number of events), cf. Fig. 1(b), and calculated the results by averaging over these ten bins. We used the same distribution also when M = B s and M = B c . For the decays with M = D s , on the other hand, we took | p Ds | = 50 GeV as a representative case, produced in LHCb by decays of B s (andB s ) mesons.
Further arithmetic averaging is made in the B-decays by averaging over the modes with B ± on the one hand and B 0 andB 0 on the other hand, because the total decay widths are somewhat different in the two cases (Γ B = 4.018×10 −13 GeV, 4.330 × 10 −13 GeV, respectively). The (partial) decay widths for the decays M → V 1 N with + 1 and − 1 are the same, as mentioned in the previous Section and explained in the Appendix. The effective detector length L in LHCb could be considered to be approximately the length of the Vertex Locator (VELO) which is about 1 m [35]. However, the length can be extended beyond the locator, to about 2.3 m [36]; we will take L = 2.3 m.
For the values of the masses and CKM matrix elements we used the central values from [33]. Similarly, we used the values of the decay constants (for the annihilation-type decays) f B = 0.1871 GeV [33], f Bc = 0.322 GeV [37] and f Ds = 0.248 GeV [37] (the average in Ref. [33] is practically the same, f Ds = 0.249 GeV). As mentioned in the previous Section, the total decay widths Γ N of N Majorana or Dirac neutrinos are obtained using the expressions of Ref. [23] (Appendix B there), based on the approach of Ref. [6].
Further, in our analysis for the rare decays with two muons in the final state we considered only the mixing elements |U µN | 2 as nonzero; for those with two taus in the final state we considered only |U τ N | 2 as nonzero. If, in addition, other mixing elements were nonzero, this would increase the total decay width Γ N , cf. Eqs. (5)-(6), and would, as a consequence, decrease the effective branching ratio Br eff of the considered rare decays and make the upper bound on the mixing element less stringent (higher).
The effective branching ratios depend on the number of produced mesons M whose rare decays we are considering. At the LHCb upgrade, the projected numbers are: N B = 4.8 × 10 12 ; N Bs = 5.76 × 10 11 ; N Bc = 2.4 × 10 10 [36]. The mesons D s are mainly produced by decays of B s mesons, with decay branching ratio around 0.9; therefore we take N Ds = 0.9N Bs = 5.18 × 10 11 .
The sensitivity upper bounds on the heavy-light mixing parameters |U µN | 2 and |U τ N | 2 at the 95% confidence limit are then obtained by requiring N events = 3.09 [45]. For example, for rare B decays we require Br eff = 3.09/(4.8×10 12 ), and obtain the corresponding upper bounds on |U µN | 2 from the absence of rare decays producing two muons, and on |U τ N | 2 when the two leptons are taus. 4 It turns out that when the two charged leptons are muons, all mentioned decays of B, B s , B c and even D s are kinematically allowed. The LNV versions of such decays are and their charge-conjugate versions. The decays of the above mesons are possible also when no vector (or pseudoscalar) mesons are produced, i.e., when the decays are of the annihilation-type. This is the case for the decays of charged When the two produced charged leptons are τ 's, the decays of the type (8) are kinematically not possible; the only decays with τ 's that are kinematically possible are The mentioned decays, Eqs. (8)- (10), are first considered to be LNV, i.e., the produced lepton pair is of equal sign (µ ± µ ± or τ ± τ ± ). In such a case, the neutrino N is considered to be Majorana. The analogous LNC decays are the decays with opposite sign lepton pairs, i.e., µ ± µ ± π ∓ gets replaced by µ ± µ ∓ π ± , and τ ± τ ± π ∓ by τ ± τ ∓ π ± . The decay rates for the analogous LNC decays are equal to those of LNV decays [34] when N is Majorana neutrino; in such a case, such LNC decays give identical sensitivity limits for the mixing parameters (|U µN | 2 ; |U τ N | 2 ) as the corresponding LNV decays. On the other hand, if N is Dirac neutrino, only LNC decays are possible, and the total decay width Γ N of Dirac neutrino is smaller than that of the corresponding Majorana neutrino (with the same mass and the same mixing parameter values). This is so because Dirac neutrinos have a significantly smaller number of decay channels than the Majorana neutrinos, resulting in about 40 per cent smaller Γ N , i.e., the coefficients N N in Eq. (6) are by about 40 per cent smaller than in the Majorana case (cf. Fig. 2 in Ref. [8]). Smaller Γ N implies that the (weakly) Γ N -dependent part of the integrand for the effective branching ratios Br eff in Eqs. (4) and (7) 1 becomes larger, and thus Br eff becomes larger. This implies that in the case of Dirac N the sensitivity limits (upper bounds) on the mixing parameters become more restrictive (smaller). This effect is expected to be appreciable only when Γ N is large, i.e., when M N is large. The resulting sensitivity limits (upper bounds) on the heavy-light mixing parameter |U µN | 2 , as a function of mass M N , for the decay processes involving muon pair are presented in Figs. 2(a)-(d). For each decay channel, we presented the limits from LNV decays (for N Majorana) and the limits when N is Dirac neutrino (LNC decays with N Dirac); the latter limits are in general only slightly better (lower) and are clearly visible only in the case of B (c) → µµπ decays at large masses M N ≈ 4.5-6.0 GeV. In these Figures we also included the present experimental upper bounds on this parameter from various experiments: the beam dump experiments PS191 [46], NuTeV [47], NA3 [48] and BEBC [49]; Belle experiment [50] which looked for the rare LNV decays B → X N (and N → π); and DELPHI experiment [51] which looked for rare ZνN decays with subsequent N decays. The present experimental LHCb upper bounds [52,53], from B meson decays with two equal sign muons, are less restrictive than the combination of Belle and DELPHI bounds.
We can see that, within the kinematically allowed ranges, the present experimental sensitivity limits can be improved the most with B → D ( * ) µµπ decays for 1.75 GeV < M N < 2.95 GeV, and by B c → µµπ decays for 2.95 GeV < M N < 5.8 GeV. Decays B s → D ( * ) s µµπ give somewhat less improved bounds than B → D ( * ) µµπ, in the mentioned mass interval [1.75, 2.95] GeV. On the other hand, the rare decays of D s , Fig. 2(d), give bounds in the region of smaller M N < 1.75 GeV (due to the smaller mass of D s ) where the present experimental bounds are more stringent (by NuTeV and PS191 experiments); the only exception is the small interval 1.75 GeV < M N < 1.85 GeV where the decay D s → µµπ gives improved upper bounds.
The results of Fig. 2(a), for the rare B-meson decays, largely agree with the results obtained by us earlier in Ref. [32], the changes appearing mainly due to the updated values of the input parameters [33]. We recall that the form factors V and A j for these decays are taken from Ref. [40] and F 1 from [39], and there they were determined in such a way that the product of them with the CKM matrix element |V cb | is independent of the value of |V cb | (this does not apply to the form factor F 0 determined in [32]). Further, here we took into account that in the B → D * µµπ decay the results are independent of the sign of charge of the µ pair; while previously in [32] we assumed a (weak) dependence on this sign. However, this difference accounts for less than one per cent in the obtained upper bounds, cf. the Appendix. Therefore, here we obtain for the sensitivity limits on |U µN | 2 from B → D * µµπ practically the same values as in Ref. [32]. From B → Dµµπ we obtain here the values lower by about 2-4 per cent than in Ref. [32], mostly due to the updated value of |V cb | in the form factor F 0 . And from B → µµπ the values of sensitivity limits |U µN | 2 are higher by about 4 per cent than in Ref. [32], mostly due to the updated value of |V ub | and of the decay constant f B .
We point out that in Ref. [32] no rare decays of B c , B s and D s were considered. Here we can see that the decays of these mesons, with two muons in the final state, can be considered as complementary to the analogous B meson decays. When two τ leptons are in the final state, the only kinematically allowed decays are B ± , B ± c → τ ± N → τ ± τ ± π ∓ . Figure 3 shows improved upper bounds on |U τ N | 2 from the decays B (c) → τ τ π. The improvement is better for B c decays, and it is in the larger kinematically allowed interval 2.0 GeV < M N < 4.4 GeV. The sensitivity limits for the case of LNC decays with N Dirac neutrino are only slightly lower, and the difference is practically invisible in the Figure. The present experimental bounds, in the presented mass range, are from DELPHI experiment [51]. The new limits are important also in view of the fact that, at the moment, the available experimental limits on |U τ N | 2 are only from one experiment (DELPHI). We point out that in Ref. [32] no such decays were considered.
Some of the considered rare decays of B c and B s appear theoretically equally or more attractive than the corresponding rare decays of B mesons. For example, the annihilation type decays B c → N → π are much less CKM-suppressed than the corresponding charged B-meson decays. Nonetheless, the obtained sensitivity limits on |U µN | 2 are only a little better in the B c case. This is so because the expected number of the produced B c mesons in LHCb upgrade is significantly lower than the expected number of the produced B mesons. Nonetheless, in the decays B, B c → τ N → τ τ π, the mesons B c give significantly better results; this is so because the pair τ τ is much heavier than µµ, and thus the fact that B c has a higher mass than B becomes important.
A similar argument applies to the decays of D s : the sensitivity limits on |U µN | 2 from D s → µµπ are less restrictive due to the suppressed expected number of produced D s .
In summary, in this work we extended the analysis of our previous work [32] where the rare decays of B mesons were considered. In the present analysis, we considered the rare decays of the mesons M = B, B c , B s and D s in the future LHCb upgrade, where in the first vertex a vector (V) or pseudoscalar (S) meson is produced together with a charged lepton 1 and an on-shell ∼ 1 GeV Majorana or Dirac neutrino N : M → V (S)W * → V (S) 1 N . The produced ∼ 1 GeV neutrino is assumed to travel and decay within the detector, N → 2 π. Further, also the annihilation-type rare decays were considered, i.e., those where no V (or S) mesons are produced, M ± → W * ± → 1 N → 1 2 π. The two charged leptons were assumed to be equal sign muons or taus: 1 2 = µ ± µ ± or τ ± τ ± , i.e., the decays were assumed to be distinctly LNV. In addition, we considered also the case of Dirac neutrino N , where the analogous decays were LNC, i.e., 1 2 = µ ± µ ∓ or τ ± τ ∓ . It turned out that, in the case that the considered decays are not detected, the upper bounds on the mixing parameters |U µN | 2 and |U τ N | 2 of ∼ 1 GeV neutrino can be significantly improved.
When the pseudoscalar M decays into a vector particle V , charged lepton 1 , and heavy (∼ 1 GeV) neutrino N , the differential decay width is significantly more complicated than when a (pseudo)scalar S is produced instead of V . Further, it is more complicated than the differential decay width when the neutrino and the charged leptons are (practically) massless. The decay M → V 1 N is presented schematically in Fig. 4. The complexity is reflected in the structure of the hadronic matrix elements; in the specific case when M = B (and V = D * ) we have and A 3 (q 2 ) denotes The first term in Eq. (A2) has a factor η = ±1, which is obtained by the application of charge conjugation operation to the hadronic matrix elements; we have η = +1 when − 1 is produced, and η = −1 when + 1 is produced. More explicitly, the (reduced) decay amplitude for the process B → D * 1 N is When making the absolute square and summing over the final leptonic helicities leads to where L µν is the lepton tensor We use for γ 5 and ε µνδη the conventions of Ref. [54]. It turns out that the η-dependence of the leptonic factor L µν and of the hadronic factor H µ H * ν cancel in the product (A5). If we define θ 1 as 6 the angle between p 1 andẑ =q , this then leads after some algebra and summation over the polarization vectors of D * to the following differential decay rate: Here, M 1 is the mass of 1 , V qQ is the relevant CKM matrix element (V qQ = V cb for M = B, B s , B c ; V qQ = V cs for M = D s ), and the following notations were used: and 7H ±1 ,H 0 andH 3 denote the following expressions containing the mentioned form factors V and A j (j = 0, 1, 2, 3): In Appendix C of Ref. [34] we already obtained the differential decay rate (A7) [Eq. (C19) there]; however, apart from some not relevant typos there, in the parametrization (A2) of the hadronic matrix elements we used the opposite sign in the first term there (proportional to V ), which is inconvenient because it then corresponds to a negative form factor V (q 2 ). As a consequence, the results of Appendix C of Ref. [34] should be reinterpreted with the substitution V → −V in the formulas there. Therefore, the term proportional to V in the differential decay rate (A7) here [i.e., the term with (H +1 ) 2 − (H −1 ) 2 ] has now the sign opposite to that in Ref. [34]. Further, in Ref. [32] we multiplied this term by η (= ±1), the structure suggested in the literature (cf.,e.g., [55,56]) which is obtained by keeping in the squared amplitude the η-dependence of the leptonic part L µν and regarding the hadronic part H µ H * ν as ηindependent. We recall the the η-dependence of the hadronic elements, Eq. (A2), is obtained by application of the charge conjugation transformation to the hadronic matrix elements. Furthermore, the measurements of the differential decay rates B + →D * 0 + 1 ν 1 [57] andB 0 → D * + − 1ν 1 [58] (for which M N ≈ 0 ≈ M 1 ) show that the differential rate is an increasing function of cos θ 1 , i.e., that the sign in front of the term ∼ cos θ 1 in Eq. (A7) is negative 8 in both types of decays, i.e., independent of η.
The use of the expression (A7) in the integration (4) is in principle sensitive to the discussed term ∼ cos θ 1 ((H +1 ) 2 − (H −1 ) 2 ), because of θ 1 -dependence of the decay probability P N there in the integrand. However, this term appears to affect the obtained upper bounds for |U N | 2 here insignificantly, by less than one per cent.
Furthermore, if we directly integrate the differential decay rate (A7) (as was the case in Ref. [34], considering P N a constant), the discussed term ∼ cos θ 1 gives exactly zero. Namely, in such a case, the integration dΩp 1 = 2πd cos θ can be performed explicitly, and the subsequent over dΩq gives then factor 4π, leading to the expression Eq. (C20b) of Ref. [34]. The explicit expression for the total decay width Γ(B → D * 1 N ) is then The explicit expression for this total decay width in Ref. [34] [Eq. (19) there] was written for the case of an approximation of the form factor A 0 as a combination of A 1 and A 2 [Eq. (17) there], but is written here, for completeness, in the form independent of this approximation.