Majorana neutrinos with effective interactions in B decays

We investigate the possible contribution of Majorana neutrinos to $B$ meson decays in an effective interaction formalism, in the mass range $0.5$ GeV $<m_N<5 $ GeV. We study the decay of the $B^-$ meson via $B^- \to \mu^- \mu^- \pi^+$ at LHCb, which is a signal for leptonic number violation and the presence of Majorana neutrinos, and put bounds on different new physics contributions, characterized by their Dirac-Lorentz structure. We also study the bounds imposed by the radiative $B$ decay ($B^- \rightarrow \mu^- \nu \gamma$) results from Belle. The obtained bounds are more restrictive than previous values found for dimension 6 four-fermion contact vectorial and scalar Majorana neutrino interactions in the context of the Left-Right symmetric model for higher Majorana masses at the LHC, showing that the direct calculation of the effective $N$ interactions contribution to different processes can help to put more stringent bounds to different UV-complete models parameterized by an effective Lagrangian.


I. INTRODUCTION.
The search for particles beyond the standard model (SM) content has been extensive in the past few years, among them sterile Majorana neutrinos N , which appear as a natural consequence in several SM extensions. The discovery of neutrino oscillations suggests that the standard neutrinos are massive particles. One of the possible ways to generate their mass is the seesaw mechanism [1][2][3][4][5], which introduces at least one right handed singlet and produces Majorana neutrinos. In this way one obtains masses for the standard neutrinos On the other hand, for smaller Yukawa couplings of order Y ∼ 10 −8 − 10 −6 , sterile neutrinos with masses around M N ∼ (1 − 1000) GeV could exist. However, in the simplest Type-I seesaw scenarios, a major drawback is that the left-right mixing parameters U lN (l = e, µ , τ ) need to be negligibly small U 2 lN ∼ m ν /M N ∼ 10 −14 − 10 −10 in order to account for light ν masses [6,7]. The mixings U lN weight the coupling of the heavy N with the SM particles, in particular with charged leptons through the V − A interaction with the active states, without making any hypothesis on the neutrino mass generation mechanism [8,9]. Such a minimal SM extension leads to contributions to LNV observables which are already close, or even in conflict, with current data from meson and tau decays, for Majorana masses M N below 10 GeV (see [8,10] and the references therein). So, also from the experimental point of view, the simple SM extensions which attribute LNV only to the mixing between heavy Majorana states and the active neutrinos are facing increasingly stringent constraints.
As suggested in [11], the detection of Majorana neutrinos (N ) would be a signal of physics beyond the minimal seesaw mechanism, and its interactions could be better described in a model independent approach based on an effective theory. One can think of an alternative treatment and consider the Majorana neutrino interactions as originating in new physics described by an unknown underlying renormalizable theory valid at a higher energy (UV) scale and parametrized at low energies by a model independent effective Lagrangian. In this approach, we consider that the sterile N interacts with the SM particles by higher dimension effective operators, taking these interactions to be dominant in comparison with the mixing with light neutrinos through the Yukawa couplings, which we neglect [12][13][14][15][16][17][18][19]. We depart from the usual viewpoint in which the mixing with the standard neutrinos is assumed to govern the N production and decay mechanisms. Here, for simplicity, we consider a scenario with only one Majorana neutrino N and negligible mixing with the ν L . among others (for comprehensive reviews see [6,20] and references therein). In particular, heavy flavor meson decays could be the place where for the first time the Majorana neutrino effects were observed or, in the absence of a discovery, this fact can be used to set limits for its coupling to SM particles. N -mediated lepton number violation in rare B meson decays has been studied, for example, in [7,8,[21][22][23][24][25][26][27][28], and the references therein. Concerning the resonant production of Majorana neutrinos in semileptonic pseudoscalar meson threebody decays, the recently measured branching ratio Br(B − → µ − µ − π + ) < 4 × 10 −9 for intermediate neutrinos with lifetimes τ N shorter than 1 ps at the LHCb experiment [29] gives the currently more stringent bounds on the mixing parameter |U µN | 2 in the case of the minimal SM extension by one Majorana neutrino (e.g. [8,30]) for Majorana masses in the range 2.5 M N 5 GeV.
In this paper we aim to exploit the recent B-decay data to constrain the possible values of the couplings that weight the contribution of different effective operators to the Majoranamediated same sign dilepton B-decay B − → µ − µ − π + and the radiative leptonic muon-mode The LHCb collaboration has presented model independent upper limits on the branching ratio of the first process [29,31], and the Belle collaboration has set new limits on the integrated differential width of the B − → µ − νγ decay [32]. The obtained bounds (for 0.5 m N 5 GeV) are more restrictive than previous values obtained for dimension 6 four-fermion contact vectorial and scalar Majorana neutrino interactions in the context of the Left-Right symmetric model for higher Majorana masses [33], and constrain the perspectives of discovery of Majorana neutrinos with effective interactions with GeV-scale masses by direct production in colliders and meson decays [17,34,35].
The paper is organized as follows. In Sect. I A we introduce the effective Lagrangian formalism. In Sect. II we calculate the B − → π + µ − µ − branching ratio and the B − → µ − νγ decay in this formalism. In Sect. III we present our results for the obtained bounds, and in Sect.IV we make our final comments.
A. Majorana neutrino effective interactions.
An appropriate way to include the Majorana neutrino into the theory is to extend the SM Lagrangian. In this work we consider an effective Lagrangian in which we include only one relatively light right handed Majorana neutrino N as an observable degree of freedom.
The new physics effects are parameterized by a set of effective operators O J constructed with the SM and the Majorana neutrino fields and satisfying the SU (2) L ⊗ U (1) Y gauge symmetry [36].
The effect of these operators is suppressed by inverse powers of the new physics scale Λ.
The total Lagrangian is organized as follows: where n is the mass dimension of the operator O (n) J . Note that we do not include the Type-I seesaw Lagrangian terms giving the Majorana and Yukawa terms for the sterile neutrinos. The dominating effects come from the lower dimension operators that can be generated at tree level in the unknown underlying renormalizable theory.
The dimension 5 operators were studied in detail in [37]. These include the Weinberg gives Majorana masses and couplings of the heavy neutrinos to the Higgs (its LHC phenomenology has been studied in [39,40]), and the operator O In the following, as the dimension 5 operators do not contribute to the studied processes -discarding the heavy-light neutrino mixings-we will only consider the contributions of the dimension 6 operators, following the treatment presented in [11]. We start with a rather general effective Lagrangian density for the interaction of right-handed Majorana neutrinos N including dimension 6 operators.
The first operators subset includes those with scalar and vector bosons (SVB), and a second subset includes the baryon-number conserving four-fermion contact terms (4-f): where l i , u i , d i and L i , Q i denote, for the family labeled i, the right handed SU (2) singlet and the left-handed SU (2) doublets, respectively. The field φ is the scalar doublet. Also γ µ are the Dirac matrices, and = iσ 2 is the antisymmetric symbol.
One can also consider operators generated at one-loop (1-loop) order in the underlying full theory, whose coefficients are naturally suppressed by a factor 1/16π 2 [11,41]: Here B µν and W I µν represent the U (1) Y and SU (2) L field strengths respectively, and σ µν is the Dirac tensor.
The effective operators above can be classified by their Dirac-Lorentz structure into scalar, vectorial and tensorial.
In this paper we will consider the B decays B − → µ − µ − π + in Sect. II B and B − → µ − νγ in Sect. II C, mediated by an on-shell Majorana neutrino N . We can thus take into account the following effective Lagrangian terms involved in the B − → µ − N and N → µ − π + processes (from eqs. (3) and (4)): The couplings α O are associated to specific operators: For the N → νγ decay, the considered Lagrangian terms come from one-loop level generated tensorial operators: , with h being the Higgs field and v its vacuum expectation value, we can write the Lagrangian (6) terms involved in our calculation (and its charge conjugate), as Here the quark fields are flavor eigenstates with family i = 1, 2, 3. In order to find the contribution of the effective Lagrangian to the B − → µ − N and N → µ − π + decays, we must write it in terms of the massive quark fields. Thus, we consider that the contribution of the dimension 6 effective operators to the Yukawa Lagrangian are suppressed by the new physics scale with a factor 1 Λ 2 , and neglect them, so that the matrices that diagonalize the quark mass matrices are the same as in the pure SM.
Writing with a prime symbol the flavor fields, we take the matrices U R , U L , D R and D L to diagonalize the SM quark mass matrix in the Yukawa Lagrangian. Thus the left-and righthanded quark flavor fields (subscript i) are written in terms of the massive fields (subscript β) as: With this notation, the SM V CKM mixing matrix corresponds to the term We first consider the lepton number violating B − → µ − µ − π + decay shown in Fig We calculate the decay of the B − meson for the two studied processes in two stages.
Firstly we obtain the decay of B − to a muon µ − and a Majorana neutrino N . Secondly, we calculate the decay of N to another muon µ − and a pion π + and the radiative decay The decay width of the B − meson is obtained in both cases in the following way with and with The details of the calculation of the total N decay width are described in [15,16]. In Figs.3a and 3b we present the results for the total width Γ N for the different sets of effective couplings, as will be described in Sect. III A.
We start with the calculation of the B − → µ − N decay width in (11). The Lagrangian massive quark fields as: where we have written the muon and massive up and b quark fields. From this Lagrangian we find the amplitude The first term in the amplitude corresponds to the W -mediated diagram which includes a SM vertex, giving the CKM V ub contribution. In the last term, we need to rearrange the field operators in order to put together the quark fields in a sandwich and the lepton fields in another. So we make a Fierz transformation taking into account a minus sign from the permutation of fermions, and then we replace it by The calculation of the leptonic matrix elements is straightforward, In order to calculate the hadronic matrix elements, we have to rely on the symmetries [42,43]. The matrix element 0|ūγ ν γ 5 b| B is a Lorentz 4-vector because the B meson is a pseudoscalar andūγ ν γ 5 b is a pseudo 4-vector. The meson state is described solely by its four momentum q µ , since it has zero spin. Therefore, q µ is the only 4-moment on which the matrix element depends and it must be proportional to q µ . Thus, we can write On the other hand, for the same reason, the matrix elements of the 4-vector, the tensor and pseudo-tensor are zero In the case of the matrix element of the scalar or pseudo-scalar interactions, we have to use the Dirac equations of motion, and we obtain the relations for the current matrix elements Putting it all together and integrating over the 2-body phase space, we obtain with M in (15) giving where The effective couplings in A V,S , as the subscript indicates, correspond to vectorial and scalar interactions.
The result is where Let us now calculate the decay N → π + µ − . According to the Lagrangian (6), the amplitude for this process can be written as The last term in (25) also needs to be modified by means of a Fierz transformation. After some algebra, it is written as In order to calculate the various factors in (25), we make use of the definition for the pion form factor and from this equation we obtain the following expressions The contribution of the pseudo-scalar quark current to the matrix element of the ordinary pion decay (27) is enhanced in comparison with the standard chirality suppressed V − A contribution and it is expected to be severely constrained by the experimental data. We finally have for the squared amplitude where 1 from which we obtain the decay width for N → π + µ − where Finally the decay width for the B − meson Γ(B − → µ − µ − π + ) is calculated according to (11) and (12), allowing us to obtain the effective branching ratio: C. The B − → µ − νγ decay.
The SM radiative leptonic B decays have been extensively studied in the literature [44][45][46][47][48][49], as they are a means of probing the strong and weak SM interactions in a heavy meson system. The measurement of pure leptonic B decays is very difficult due to helicity suppression and the fact of having only one detected final state particle. On the other hand the radiative modes, with an extra real final photon, can be even larger than the pure leptonic modes as they escape helicity suppression and are also easier to reconstruct.
The Belle collaboration has recently released an analysis of the full Belle experiment dataset [32] using new theoretical inputs [49] for the QCD calculations and new algorithms prepared for the Belle II experiment. They obtain the experimental bound ∆Br exp B − →µ − νγ < 3.4 × 10 −6 for the integrated partial branching ratio of the muon-mode radiative B decay.
We consider the SM and the effective contribution coming from the B → µN followed by N → νγ reaction, and use the Belle bound to set limits on the one-loop generated effective couplings involved in this last decay mode, as will be discussed in Sect. III A.
The SM B → µνγ differential decay width in the B meson rest frame can be parameterized as [48]  We call ∆Br SM to the integrated partial branching ratio in the energy range Here for kinematic reasons E max γ = m B /2 and the minimal photon energy infrared cutoff E cut is such that the theoretical QCD treatment remains valid. As we will use the 2 These values are also consistent with the central values given in figures (7) and (8) of reference [49], for the inverse moment of the leading twist light cone distribution amplitude λ B value given by Belle [32].
latest Belle results for the experimental limit ∆Br exp , we take E cut = 1 GeV, as in ref. [32]. The value we obtain for our estimation of the partial branching ratio in the SM is ∆Br SM ∼ 5 × 10 −7 , which is of the order of the values recently considered in ref. [50].
Here M B→µN is the amplitude presented in (22) and M N →νγ is the amplitude of the radiative N → νγ decay allowed by the one-loop generated operators (5) in the Lagrangian (8), again corresponding to the second fermion family i = 2: Thus, multiplying and dividing (34) by the partial width Γ (N →νγ) we have where Br(N → νγ) is the branching ratio in (13). Partially integrating the phase space, the last factor in (36) can be written as where x = k 0 /m N , with k 0 the energy of the photon in the Majorana N rest frame. The distribution in the B meson rest frame is obtained by a Lorentz transformation. Here, as in (33) E γ is the photon energy in the B rest frame, so We use (38) in order to transform the distribution Thus, for −1 < cos θ < 1 we have 1 Integrating in x and cos θ we have In order to obtain the partial branching fraction for B → µνγ with photon energy E γ > E cut we integrate (40) The integration region is shown in Fig. 5 for E cut = 1 GeV.

A. Numerical treatment
The aim of this work is to study the bounds that can be set on the different couplings α J in the effective dimension 6 Lagrangian (2) involved in N mediated B decays by exploiting the experimental results existing on the B − → µ − µ − π + [31] and B − → µ − νγ [32] processes.
which is derived from the contribution of the operator O (3) and allows a direct comparison with the mixing angles in the Type I seesaw scenarios [7].
Some of the operators involving the first fermion family (with indices i = 1) are strongly constrained by the neutrinoless double beta decay bounds, currently obtained by the KamLAND-Zen collaboration [51]. Following the treatment already made in [16], the values of the 0νββ-decay constrained couplings α (1) The B to final muon decays studied in this work allow us to set bounds on the couplings involving the second fermion family (generically α J ). As we found in Sect. II, the B → µ − N effective decay depends on the couplings appearing in the vectorial (A V ) and scalar (A S ) interactions in (23), the N → µ − π + depends on the C V and C S couplings in (29) and the N → νγ depends on the one-loop tensorial couplings in (35). In the first two cases, new quark flavor-mixing matrices U R,L and D R,L appear.
For the numerical treatment we name the quark mixing matrix combinations in (23) as The corresponding notation for the combinations appearing in C V and C S in (29) is set by changing the superscripts ub for ud (and 3 → 1 in the D matrices superscripts), as they involve the d quark instead of the b quark.
All these new mixing matrices are unknown, and in principle their entries may be found by independent measurements, as is done in the case of the SM V CKM matrix. In this occasion we will make an ansatz and consider that all the Y ub values in (43) shall be of the order of the SM V ub value, taking it as a measure of the strength of the coupling between the u and b quarks. Correspondingly, we will consider the Y ud values to be of the order of the SM V ud CKM mixing.
This allows us to consider A V and A S in (24) and C V and C S in (30) for the numerical treatment as and set bounds on the possible values of these effective couplings using the B-decay data.
As we would like to disentangle the kind of new physics contributing to the Majorana neutrino interactions, for the numerical analysis we will consider different benchmark scenarios for the effective couplings, where we switch on/off the operators with distinct Dirac-Lorentz structure: vectorial, scalar and the tensorial one-loop generated operators. If we call (V, S, L) the factors multiplying the vectorial, scalar and one-loop generated operators respectively, we can define six sets, presented in the table I.
In order to exploit B-decay data to put bounds on the effective couplings in table I, we will take them as equal to the same value α, and use the experimental results constraining the value of the combination U 2 defined in (42). We have α = 2Λ 2 /v 2 √ U 2 for the tree-level generated operators (which are the vectorial and scalar operators), and in the case of one-loop generated operators we have α = 1 This allows us to write the numerical results for the total Majorana neutrino decay width Γ N , the branching ratio in (31) and the integrated effective branching ratio ∆B ef f r (B − → µνγ) in (41) as a function of the Majorana neutrino mass m N and the U 2 combination.
The sets 1, 2 and 3 in Tab. I take into account the contributions of the one-loop generated effective couplings in (8) to the N decay width. In particular this sets allow for the existence of the N → νγ decay channel represented in Fig. 2. As we found in [16], this channel gives the dominant contribution to the N decay width for the low mass m N range considered in this work. The sets 4, 5 and 6 discard this contribution. As can be seen in Fig. 3, the total Γ N width is around three orders of magnitude higher in the sets 1, 2 and 3 (Fig.3a) than in the sets 4, 5 and 6 ( Fig.3b). In fact, as the scalar and vectorial couplings contribution to the N decay in this mass range is so poor, the three curves in Fig.3a cannot be distinguished in the plot scale. This effect in the Γ N value will explain many of the differences in the bounds we obtain for the U 2 combination when we consider one group of sets or the other, as will be discussed below.

B. Obtained bounds
We start by discussing the bounds obtained from the LHCb results on the B − → µ − µ − π + decays. The LHCb collaboration has presented a search for Majorana-mediated B − → µ − µ − π + decays [29], where they obtain model independent limits on the branching ratio Following the procedure in ref. [52], we convert the model independent LHCb upper limits on the branching ratio Br(B − → µ − µ − π + ) into limits on the combination U 2 defined in (42). For each value of m N (which fixes the value of τ N for a given U 2 value in the effective model) we scan through the values of U 2 for which our computed branching fraction (31) equals the upper bound in ref. [29]. For the experimental values, we consider the data in Fig.   6, which reproduces the values presented in Figure 5 of ref. [29]. The obtained constraints on the U 2 values are presented in Fig. 7a for the coupling sets 1, 2 and 3 in Tab. I, and As can be seen in the plots, the bounds we obtain for U 2 in the sets 1, 2 and 3 ( Fig.   7a) are weaker than those we get in the sets 4, 5 and 6 ( Fig. 7b). This is explained by the different Γ N values in the two groups of sets discussed above: in the case of the sets 4, 5 and 6, the value of the branching ratio Br ef f (B − → µ − µ − π + ) in (31) is around 10 3 times higher, because the Γ N factor in the denominator is lower than in sets 1, 2 and 3, and thus we obtain more restrictive bounds for the U 2 combination when we do not take into account the one-loop generated effective couplings contribution.
On the other hand, among the sets in each figure, we find that we place stronger bounds on the scalar couplings (considering their sole effect in set 6, and with one-loop couplings in set 3). This is due to the presence of the light quark masses in the denominators of the scalar terms, enhancing these contributions to the Br(B − → µ − N ) in (24) and the  [29], considering the effective coupling sets defined in Table I. The black full line curve represents the revised bound presented in [52] for |U µN | 2 .  (41), with the Belle result [32].
We now scan for the values of U 2 for which the complete theoretical value ∆B r = ∆B SM r + ∆B ef f r equals the upper bound ∆Br exp B − →µ − νγ < 3.4 × 10 −6 , for each mass m N . The bounds we obtain for U 2 from this procedure are presented in Fig. 8.
As the one-loop generated operators need to be non-zero to allow the B − → µ − νγ decay, bounds are established just for the sets 1, 2 and 3. Again, we obtain stronger bounds on the scalar operators, due to their contribution to the N production in B decay (24). These bounds are compatible with and less restrictive than the ones obtained from the LNV process B − → µ − µ − π + in Fig. 7a. When one-loop generated interactions are also taken into account, the bound is relaxed to α ≤ 0.26 in sets 1 and 3.
The α ≤ 0.16 bound should be compared for instance with the upper bound our group considered for the calculation of the contribution of scalar and vectorial effective Majorana interactions to the LNV same-sign dilepton signal pp → µ + µ + jj in the LHC [17]. In that early work we estimated an upper bound α ≤ 0.3 coming from the heavy neutrinos search results at Belle [53]. Other works also took into account the same bound for the calculation of prospects for the observation of e + e − → νN → νγ at Belle-II and the ILC [34]. The revision of these results is left for future work.   The result is obtained performing a reinterpretation in terms of the LRSM model of the LHC limits from heavy Majorana neutrino direct production at √ s = 8 TeV in the dilepton channel pp → W * R → N µ ± → µ ± µ ± + nj [54]. The most stringent bounds on α V 0 ,S 3 are obtained considering Br(N → µX) ∼ 1, so that the N decays preferably to muons. These are (taking Λ = 1 TeV for comparison) α V 0 0.23 and α S 3 0.45 for m N = 100 GeV.
The comparison suggests that the direct calculation of the effective N interactions contribution to different processes can help to put more stringent bounds to different UV-complete models parameterized by the effective Lagrangian in (2).

IV. FINAL REMARKS
We have considered heavy Majorana neutrinos coupled to the ordinary matter in a general way by dimension 6 effective operators satisfying the SM electroweak symmetry. According to these interactions these neutrinos would be produced in the decay of B mesons, and subsequently decay to standard particles. In particular, we exploited the non-observation of the B − → µ − µ − π + decay in the LHCb [29] and put limits to the couplings of the different effective operators contributing to this decay in the Majorana mass range 0.5 GeV < m N < 5 GeV. These upper bounds are presented in Fig. 7. Also for this m N range, we have considered the bounds coming from the radiative decay B − → µ − νγ by the Belle experiment [32]. This allows us to set bounds directly on the one-loop generated operators. These bounds are compatible with and weaker than the ones we derive form the LNV process B − → µ − µ − π + and are shown in Fig. 8.
The obtained bounds (for 0.5 m N 5 GeV) are more restrictive than previous values obtained for dimension 6 four-fermion contact vectorial and scalar Majorana neutrino interactions in the context of the Left-Right symmetric model for higher Majorana masses [33].
The comparison suggests that the direct calculation of the effective N interactions contribution to different processes can help to put more stringent bounds to different UV-complete models parameterized by the effective interaction formalism. The obtained upper bounds also constrain the perspective of discovery of Majorana neutrinos with GeV-scale masses by direct production in colliders and meson decays [17,34,35].