Lepton number violation in D meson decay

The lepton number violating process can be induced by introducing a fourth generation heavy Majorana neutrino, which is coupled to the charged leptons of Standard Model. There have been many previous studies on the leptonic number violating decay processes with this mechanism, we follow the trend to study the process: $D \to K \ell \ell \pi$ with the same-sign dilepton final states. We restrict ourself to certain neutrino mass regions, in which the heavy neutrino could be on shell and the dominant contribution to the branching fraction comes from the resonance enhanced effect. Applying the Narrow Width Approximation, we found that upper limit for the branching fractions for $D^0 \to K^-\ell^+\ell^+\pi^-$ are generally at the order of $10^{-10}$ to $10^{-9}$, if we take the most stringent upper limit bound currently available in the literature for the mixing matrix elements. We also provide the constrains, which is competitive compared to the LNV B decays, on the mixing matrix element $|V_{eN}|^2$ based on the upper limit of $D^0 \to K^- e^+ e^+ \pi^-$ estimated from Monte-Carlo study at BESIII.


I. INTRODUCTION
The discovery of neutrino oscillations [1][2][3] and observation of unexpected large θ 13 [4] have convincingly shown that neutrinos have finite mass and that lepton flavor is violated in neutrino propagation. The generation of neutrino masses is still one of the fundamental puzzle in particle physics. To obtain the non-vanishing mass of neutrinos, one must make a minimal extension of the Standard Model (SM) by including three right-handed neutrinos.
Once the consensus that the neutrino is a massive fermion have been reached, another urgent task is to figure out whether the neutrinos are Dirac or Majorana particles, the latter case is characterized by being their own antiparticles.
The fact is, at present, the Majorana neutrino is one of the favorite choices for the most of theories, since the masses of the observed light neutrinos could be naturally derived from heavy neutrinos via the so-called "see-saw" mechanism [5][6][7][8][9][10][11]. Owing to the new heavy neutrinos Majorana nature, it is its own antiparticle, which allows processes that violate lepton-number conservation by two units. Consequently, searches for Majorana neutrinos are of fundamental interest.
There is one promising method to probe the Majorana nature of neutrinos, i.e., neutrinoless double beta (0νββ) decay [12][13][14][15][16][17][18][19][20][21][22][23][24]. As we have mentioned above, the existence of Majorana mass term could induce the lepton number violation (LNV) by exchange of virtual Majorana neutrinos between two associated beta decays. Although the first double-beta decay was proposed as early as 1935 by Goepper-Mayer, it was until four years later Furry first calculated the (0νββ) decay based on the Majorana theory [25]. These early exploration gave an impetus to many years of experimental and theoretical research. It is an interesting question that whether the moderately heavy sterile neutrino with a mass from few hundred MeV to a few GeV exists. If such neutrino exists, the decay rates of these processes can be substantially enhanced by the effect of neutrino-resonance which is induced by taking the mass of virtual Majorana neutrino as a sterile neutrino mass. As a consequence, we may have an chance to observe such LNV processes, or set up stringent upper bounds on sterile neutrino mass and mixing matrix elements. Thanks to the important significance, these LNV processes have been extensively studied in both theory [26] and experiment [27][28][29][30][31][32][33][34].
In recent years, most of the studies of these LNV processes focus on the three-body and four-body ∆L = 2 decays of K, D and B mesons, as well as tau lepton decays. The decay processes: K + , D + , B + → ℓ + ℓ + M − , where M denotes a vector or pseudoscalar meson and ℓ = e, µ, τ , have been extensively studied in [26,[35][36][37][38][39]. In the meanwhile, the tau lepton decay τ − → ℓ + M − 1 M − 2 has also been discussed in [37,[40][41][42][43]. Compared to the three-body LNV processes, the studies of four-body LNV decays are relatively rare. The B → Dℓℓπ four-body decay has been recently calculated in [44], and the four-body decays [45]. Recently, more and more attentions are directed to the four-body LNV decay since these processes can give more abundant kinematical regions than the three-body LNV decays, thus we can grope more wide range of the resonance neutrino mass. In this paper, we will study the four-body LNV decays D → Kπℓ + ℓ + .
In Particle Data Book (PDG) [34], the upper limits for the branching fractions of threebody D → ℓ + ℓ + M − decays are at the level of 10 −6 to 10 −4 . For the four-body LNV D decays, we have less information until recently. In [46], the branching fractions of D meson decay D 0 → M − 1 M − 2 l + l + were reported to be the order 10 −5 ∼ 10 −4 . From experimental side of view, the BESIII experiment is taking data at open-charm threshold, and an integrated luminosity of 2.9 fb −1 data sample has been collected at ψ(3770) by the BESIII detector [47,48]. Estimated from Monte-Carlo (MC) sample of the same luminosity, the sensitivity of D → Kπℓ + ℓ + can reach the level of 10 −9 [49]. It is interesting to investigate D → Kπℓ + ℓ + decays mediated by an on-shell Majorana neutrino.
The paper is organized as follows, in Section II, we give a brief introduction to the theoretical framework involving the LNV decays; in Section III, we sketch some techniques in our computations, like phase space parametrization and Narrow Width Approximation, etc.; in Section IV, we give our numerical results for the branching fractions of LNV D meson decays; and finally comes the summary, we also provide the analytical results of the squared amplitude in the Appendix A.

II. THEORETICAL FRAMEWORK
We consider the LNV four-body decay of D meson: where D with momentum p and K with momentum p 1 can be charged or neutral, two ℓ + are charged leptons with momenta p 2 and p 3 , respectively, and the charged pion π − has dd c d Following the previous studies [26,50], such LNV process can be induced through a Majorana neutrino N coupled to the charged leptons ℓ, such gauge interactions are described by the following vertex in the Lagrangian: where P L = 1 2 (1 − γ 5 ), N is the mass eigenstate of the fourth generation Majorana neutrino, V ℓN is the mixing matrix between the charged lepton ℓ neutrino ν ℓ and heavy Majorana neutrino N, their restrictive bounds are reported in [51].
At leading order, the Feynman diagrams are displayed in Fig. 1 for the charged and neutral D meson decays, the diagrams with the charged leptons exchanged are also included.
If the neutrino mass is from a few hundred MeV to several GeV, the neutrino propagator in diagrams (a), (c), (e) and (f) in Fig. 1 could be on-shell, i.e., the neutrino becomes a resonance, and its contribution to the decay can be much enhanced due to such neutrinoresonance effect, so the contributions from other diagrams are negligible compared to such enhanced diagrams, since the neutrino can not become a resonance in those diagrams due to the kinematic restrictions. In addition, we note that diagrams (c) and (f) are suppressed The transition amplitude for such LNV process can be written as: GeV −2 is the Fermi constant, and V cs and V ud are the CKM matrix element. We have already replaced the neutrino propagator with its resonant type: and q N is the momentum of heavy Majorana neutrino, q ′ N is the same except with two charged leptons exchanged, m N is the mass of such heavy Majorana neutrino, Γ N is the total decay width of the heavy Majorana neutrino.
The decay width of the heavy Majorana neutrino can be obtained by adding up all contributions of neutrino decay channels which can be opened up at the mass m N [26]: where m s in argument of Heaviside θ function are the masses of final state particles in the corresponding decay channel. All expressions for these decay width can be found in Appendix C of [26].
The matrix element involving π in Eq.(4) is related to the decay constant of charged π by: where f π is the decay constant of the charged π. Then squared amplitude can be obtained as: where the factor PN 1 and PN 2 are defined as: The complete analytical expression for the squared amplitude can be found in Appendix A.

A. Kinematics for four-body decay
The kinematics of the four-body decay D → Kℓ + ℓ + π − can be described in terms of five independent variables: {s 12 , s 34 , θ 12 , θ 34 , φ}, which have the following geometrical definitions [52] (see Fig. 2): (i) s 12 ≡ (p 1 + p 2 ) 2 is the square of invariant-mass of Kℓ system; (ii) s 34 ≡ (p 3 + p 4 ) 2 is the square of invariant-mass of πℓ system; (iii) θ 12 is the angle between the K three-momentum in the Kℓ rest frame and the line of flight of the Kℓ in the D rest frame; (iv) θ 34 is the angle between the π three-momentum in the πℓ rest frame and the line of flight of the πℓ in the D rest frame; (v) φ is the angle between the normals to the planes defined in the D rest frame by the Kℓ pair and the πℓ pair.
The angular variables are shown in Fig. 2, where K is the K three-momentum in the Kℓ center-of-mass(CM) frame and π is the three-momentum of the π in the πℓ CM frame. Let v be the unit vector along the Kℓ direction in the D rest frame,ĉ the unit vector along the projection of K perpendicular tov, andd the unit vector along the projection of π perpendicular tov. We have , cos φ ≡ĉ ·d, sin φ ≡ (ĉ ×v) ·d .
The Lorentz invariant phase space for the four-body decay is defined as where q ij = p i + p j andβ andβ ij are defined as: Here we decompose a four-body phase space integral into a product of two-body phase space integrals. This is very useful if one considers a production of two particles, each of which subsequently decays into two-body state.

B. Narrow Width Approximation
With the squared amplitude in Eq.(8) and phase space parametrization in Eq.(11), we are ready to get the decay rate for D → Kℓℓπ using the decay rate formula: Let's look at the Eq.(4) again before we start to perform the phase space integral. In general, q N = q ′ N , and it is convenient to split up the individual resonant contributions by the Single-Diagram-Enhanced (SDE) multi-channel integration method [53]. To do this, we define the where each M i corresponds to the amplitude for a single diagram, then the total amplitude squared is given by The amplitude squared now splits up into the functions f i defined above and the phase space integration can be done for each f i separately, moreover the peak structure of each f i is the same as of the single squared amplitude |M i | 2 . When the width Γ N of the heavy Majorana neutrino is very small compared to the neutrino mass m N , we can apply the Narrow Width Approximation(NAW): Applying NAW and SDE multi-channel integration method, we can make convenient simplification for the phase space integration and the computation can be carried out in parallel.
The contribution from each f i can be added up after phase space integration.

IV. NUMERICAL RESULTS
The matrix element K|sγ µ (1 − γ 5 )c|D for D to K can be parameterized as: where q = p − p 1 , and f + and f 0 are two form factors. When we take the lepton mass as zero, the f 0 (q 2 ) will not contribute, the strong interaction dynamics can be described by a single form factor f + (q 2 ). We use the Modified Pole (MP) [54] ansatz to parameterize the the form factor f + (q 2 ): where m pole is the pole mass which is predicted to be the D * − mass, and α pole is a free parameter. We take the parameters from CLEO-c measurement [55] as follows:  As for the neutrino mixing matrix elements, we take the most stringent upper limit currently available in the literature [26,43] in our numerical computation, these values are listed in TABLE III for both |V eN | 2 and |V µN | 2 .  .3 and FIG.4 in [26,43]. Since |V τ N | 2 is not controlled by experimental limits due to the absence of experimental data on the LNV or LFV processes involing two τ -leptons, we use an ad hoc assumption frequently used in the literature in our computation for total decay width of heavy majorana neutrino Γ N . With this assumption, it is easy to derive the following relation between decay rates Γ of D → Kπℓℓ with different mixing matrix elements according to the NAW approximation in Eq. (16): To make the heavy mojorana neutrino to be resonant, the neutrino mass m N should satisfy the following kinematical restrictions: We take the neutrino mass m N from 150 MeV to 1000 MeV, the numerical results for the upper limits of branching fractions are listed in the TABLE IV and V. We find that the The masses of heavy Majorana neutrino are in unit of MeV, and the total decay width of D 0 is Γ tot = 1.605 × 10 −9 MeV. magnitude of the branching fraction for the D 0 → K − e + e + π − and D 0 → K − µ + µ + π − are TABLE V: Upper limits on the branching fractions for D + →K 0 + ℓ + + ℓ + + π − and D + → K 0 + µ + + µ + + π − . The masses of heavy Majorana neutrino are in unit of MeV, and the total decay width of D + is Γ tot = 6.329 × 10 −10 MeV. At the BESIII experiment, estimated from 2.9 fb −1 MC sample, the upper limit for the decay of D 0 → K − e + e + π − obtained in Ref. [49] is about 1.0 × 10 −9 at 90% confidence level. It is already below the upper limit with the most stringent mixing matrix elements available in the literature for the low majorana neutrino mass m N ≤ 500MeV. In Fig. 3 we plot the exclusion regions provided by this MC study. From such estimation, we expect the BESIII experiment can provide competitive constraints on the mixing matrix element |V eN | 2 compared to the LNV B decays [33]. In the future, at charm factory, about 1 ab −1 integrated luminosity will be collected per year [56], and an improvements by three orders of magnitude on the branching fractions would yield more stringent constrain on both the mass of Majorana neutrino and the mixing matrix elements.

V. SUMMARY
We have investigated the LNF four-body decay of D meson by introducing a fourth generation heavy Majorana neutrino N, which is coupled to the charged leptons from the SM, the branching fractions for such decays are dependent on the heavy neutrino mass m N , we found that upper limit of the branching fractions are generally at the order of 10 −10 to 10 −9 , if we take the most stringent limit bound currently available in the literature for the mixing matrix element. We also provide competitive constrains on the mixing matrix element |V eN | 2 based on the upper limit of D 0 → K − e + e + π − estimated from MC study at BESIII.
Note added: After the calculation was completed and while we were preparing the draft, a related work recently appeared in arXiv [57], which also investigated the leptonnumber-violating D meson four-body decay processes. Aside from the different strategy in parameterizing the D → K form factor (The authors of [57] used the Bethe-Salpeter approach to estimate those form factors, while we use the Modified Pole ansatz whose parameters are directly extracted from the latest published CLEO-c data [55]), our numerical predictions for the branching ratios are in agreement with theirs in magnitude.