Search for heavy Majorana neutrino in lepton number violating decays of D->K pi e+ e+

Using the data sample of an integrated luminosity of 2.93 fb$^{-1}$ taken at the center-of-mass energy of 3.773 GeV, we search for the Majorana neutrino in the lepton number violating decays $D\to K \pi e^+ e^+$. No significant signal is observed, and the upper limits on the branching fraction at the 90\% confidence level are set to be $\mathcal{B}\,(D^0 \to K^- \pi^- e^+ e^+)<2.7\times10^{-6}$, $\mathcal{B}\,(D^+ \to K_S^0 \pi^- e^+ e^+)<3.3\times10^{-6}$ and $\mathcal{B}\,(D^+ \to K^- \pi^0 e^+ e^+)<8.5\times10^{-6}$. The Majorana neutrino is searched for with different mass assumptions ranging from 0.25 to 1.0 GeV/$c^2$ in the decays $D^0\to K^- e^+ \nu_N(\pi^- e^+)$ and $D^+\to K_S^0 e^+ \nu_N(\pi^- e^+)$, and the upper limits on the branching fraction at the 90\% confidence level are extracted to be at the level of $10^{-7} \sim 10^{-6}$, depending on the mass of Majorana neutrino. The constraints on the mixing matrix element $|V_{eN}|^2$ are also evaluated.

Using the data sample of an integrated luminosity of 2.93 fb −1 taken at the center-of-mass energy of 3.773 GeV, we search for the Majorana neutrino in the lepton number violating decays D → Kπe + e + . No significant signal is observed, and the upper limits on the branching fraction at the 90% confidence level are set to be B (D 0 → K − π − e + e + ) < 2.7×10 −6 , B (D + → K 0 S π − e + e + ) < 3.3×10 −6 and B (D + → K − π 0 e + e + ) < 8.5 × 10 −6 . The Majorana neutrino is searched for with different mass assumptions ranging from 0.25 to 1.0 GeV/c 2 in the decays D 0 → K − e + νN (π − e + ) and D + → K 0 S e + νN (π − e + ), and the upper limits on the branching fraction at the 90% confidence level are extracted to be at the level of 10 −7 ∼ 10 −6 , depending on the mass of Majorana neutrino. The constraints on the mixing matrix element |VeN | 2 are also evaluated.

I. INTRODUCTION
In the Standard Model (SM), due to the absence of right-handed neutrino component and requirements of SU(2) L gauge invariance and renormalizability, neutrinos are postulated to be massless. However, the observations of neutrino oscillation phenomena [1][2][3][4] have convincingly shown that neutrinos are of a very tiny mass, which provides the first evidence for physics beyond the SM. Theoretically, the leading model to accommodate the neutrino masses is the so-called "see-saw" mechanism [5][6][7][8], in which the mass (m ν ) of neutrino can be formulated by m ν ∼ y 2 ν υ 2 /m N , where y ν is a Yukawa coupling of neutrino to the Higgs field, υ is the Higgs vacuum expectation value in the SM, and m N is the mass of a new massive neutrino state N . The smallness of m ν can be attributed to the existence of the new neutrino state N with high mass.
In the "see-saw" mechanism, the SM neutrinos turn out to be Majorana particles, who are of its own antiparticle. The effects of Majorana neutrino can be manifest through the process violating the lepton-number (L) conservation by two units (∆L = 2). Consequently, experimentally searches for the Majorana neutrinos manifested through the process of lepton-number violating (LNV) ∆L = 2 are of great interest. Different ∆L = 2 processes at low and high energies have been proposed in the literature [9]. Among them, an interesting source of LNV processes is given by exchanging a single Majorana neutrino with a mass of the order of heavy flavor mass scale [9]. The effects of such a heavy neutrino with the mass around 100 MeV/c 2 to a few GeV/c 2 have been widely studied in ∆L = 2 three-body and four-body of heavy flavor mesons or τ lepton decays by different experiments, as summarized in Ref. [10], but no evidence has been seen so far. The ∆L = 2 processes of D mesons have been reported by E791 collaboration [11] with the upper limits (ULs) of the decay branching fraction (BF) ranging 10 −5 ∼ 10 −4 .
In this paper, we present the first studies of LNV process with ∆L = 2 in D meson decays D 0 → K − π − e + e + , D + → K 0 S π − e + e + and D + → K − π 0 e + e + . These precesses can occur through Cabibbo-favored (CF) and doubly Cabibbo-suppressed (DCS) decays by mediating a Majorana neutrino [12], as depicted in Fig. 1. As usually, the DCS processes (Figs. 1(c) and 1(d)) are expected to be suppressed by a factor |V cd V us /V cs V ud | ∼ 0.05 [13] with respect to the CF processes (Figs. 1(a) and 1(b)). In this analysis, we search for the LNV process with ∆L = 2 for the above three processes as well as the Majorana neutrino with different m N hypotheses in the CF processes. Additionally, the constraints on the mixing matrix element |V eN | 2 are also estimated m N dependently. The analysis is carried out based on the data sample with an integrated luminosity of 2.93 fb −1 at the center-ofmass (C.M.) energy ( √ s) of 3.773 GeV collected with the BESIII detector. Throughout the text, the charged conjugated modes are always implied implicitly.

II. DETECTOR AND MONTE CARLO SIMULATION
The BESIII detector is a magnetic spectrometer [14] located at the Beijing Electron Positron Collider (BEPCII) [15]. The cylindrical core of the BESIII detector consists of a helium-based multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over 4π solid angle. The chargedparticle momentum resolution at 1 GeV/c is 0.5%, and the dE/dx resolution is 6% for the electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end cap) region. The time resolution of the TOF barrel part is 68 ps, while that of the end cap part is 110 ps. Simulated samples produced with the geant4based [16] Monte Carlo (MC) package which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the backgrounds. The simulation includes the beam energy spread and initial state radiation (ISR) in the e + e − annihilations modelled with the generator kkmc [17].
The cocktail MC sample consists of the production of DD pairs with consideration of quantum coherence for all neutral D modes, the non-DD decays of the ψ(3770), the ISR production of the J/ψ and ψ(3686) states, and the continuum processes incorporated in kkmc [17]. The known decay modes are modelled with evtgen [18] using BFs taken from the Particle Data Group [13], and the remaining unknown decays from the charmonium states with lundcharm [19]. The final state radiations (FSR) from charged final state particles are incorporated with the photos package [20]. The cocktail MC sample is generated largely to study the possible background sources, and normalized to the luminosity of the data sample in the analysis.
To study the detector efficiencies of the LNV ∆L = 2 processes, the signal D meson is assumed to decay uniformly in phase space, while in searching for the Majorana neutrino, the exclusive MC samples D 0 → K − e + ν N and D + → K 0 S e + ν N with ν N → π − e + are generated with different m N assumption, where the angular distributions are simulated according to the squared amplitude in Eq. (8) of Ref. [12]. The form factor is described with modified pole approximation.

III. EVENT SELECTION
Charged tracks in a candidate event are reconstructed from hits in the MDC. The charged tracks other than those from K 0 S decay are required to be within 10 cm of the interaction point (IP) in the beam direction and within 1 cm in the plane perpendicular to the beam, as well as satisfying | cos θ| < 0.93, where θ is the polar angle relative to the beam direction. The information of TOF and dE/dx is combined to evaluate particle identification (PID) probabilities (prob) for the π and K hypotheses, a π (K) is identified by requiring prob(π) > prob(K) (prob(K) > prob(π)). To identify electron or positron, besides TOF and dE/dx, the information from the EMC is included additionally to compute probability. An electron or positron is required to satisfy prob(e)/(prob(e) + prob(π) + prob(K)) > 0.8, and E/p > 0.8, where E and p are the deposited energy in the EMC and the track momentum measured in the MDC, respectively.
The K 0 S candidates are reconstructed with a vertexconstrained fit for pairs of oppositely charged tracks, which are required to be within 20 cm of the IP along the beam direction, but without constraint in the transverse plane, where the two charged tracks are assumed to be pions without PID. A vertex fit is carried out to make sure that the two selected tracks are originated from a common vertex by requiring the corresponding fit χ 2 less than 100. The resulting decay vertex is required to be away from the IP by greater than 2 times of the resolution. The K 0 S candidates are further required to have an invariant mass within [0.487, 0.511] GeV/c 2 .
Electromagnetic showers are reconstructed from clusters of energy deposited in the EMC. The energy deposited in nearby TOF counters is included to improve the reconstruction efficiency and energy resolution. The photon candidate showers must have a minimum energy of 25 MeV in the barrel region (| cos θ| < 0.80) or 50 MeV in the endcap regions (0.86 < | cos θ| < 0.92). To suppress showers originated from charged particles, a photon must be separated by at least 10 • from any charged tracks. To suppress electronic noise and energy deposits unrelated to the event, timing information from the EMC for the photon candidates must be in coincidence with collision events i.e., 0 ≤ t ≤ 700 ns. The π 0 candidates are reconstructed from pairs of photons. Due to the worse resolution in the endcaps of the EMC, π 0 candidates reconstructed with two photons at the endcaps of the EMC are rejected. The invariant mass of two photons is required to be within [0.115, 0.150] GeV/c 2 for π 0 candidates. In the following analysis, the photon pair is kinematically constrained to the nominal mass of the π 0 to improve the resolution of π 0 momentum.
In order to improve the degraded momentum resolution of the positron due to the effects of FSR and bremsstrahlung, the photon, which is of the energy greater than 30 MeV, separated more than 20 • apart from any showers in the EMC originated from the charged tracks, and within a cone of 5 • around the positron direction before curvature in the magnetic field, is included additionally for the positron kinematic, namely FSR recovery process.
In the analysis, the signal candidates of D meson LNV decay are searched for by performing a single tag (ST) method. Two variables, the beam energy constrained mass M BC and the energy difference ∆E,  Table I. The potential background sources are examined by the cocktail MC sample. The dominant contributions are from the processes ψ(3770) → DD with D → Keν e due to large BFs and the processes e + e − → qq, but no peaking background is observed on the M BC distribution.

IV. SIGNAL EXTRACTION
The signal yields are extracted by performing the unbinned maximum likelihood fit on the M BC distribution of survived candidate events. In the fit, the background shape is described by an ARGUS function [21], and the signal shape is modeled by the MC simulated shape convolved with a Gaussian function which compensates the resolution difference between data and MC simulation. In practice, the width of Gaussian function is fixed to be 0.32 MeV/c 2 , obtained by studying the control sample D 0 → K − π + π + π − decay. The fitting curves are shown in Fig. 2. The BFs, B D→Kπe + e + , are calculated by where N sig is the signal yield extracted from the fit, N tot DD is the total number of DD pairs, which are (8, 296 ± 31 ± 65) × 10 3 for D + D − pairs and (10, 597 ± 28 ± 98) × 10 3 for D 0D0 pairs [22], ǫ is the detector efficiency, obtained from the corresponding MC simulation, and B is the decay branching fraction of intermediate state, i.e., 1 in the decay D 0 → K − π − e + e + due to no intermediate state involved, B K 0 S →π + π − in the decay D + → K 0 S π − e + e + and B π 0 →γγ in the decay D + → K − π 0 e + e + , where B K 0 S →π + π − and B π 0 →γγ are taken from the world average values [13]. A factor 2 in the denominator indicates both D andD mesons in every event are included.
Since no obvious signal is observed, the ULs at the 90% confidence level (CL) on the BFs of D → Kπe + e + decays are set after considering the effect of systematic uncertainties.

V. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties arise from the several aspects including tracking and PID efficiencies of charged tracks, K 0 S and π 0 reconstructions, the total number of DD pairs, cited BFs of K 0 S → π + π − and π 0 → γγ decays, ∆E requirement, FSR recovery, modeling for detector efficiency and fitting M BC .
Systematic uncertainties from tracking efficiency of K, π and e are assigned to be 1.0% per track [23,24]. For the PID efficiency, the systematic uncertainties for K(π) and e are 0.5% and 1.0% per track [23,24], respectively. Systematic uncertainties from K 0 S and π 0 reconstruction are taken to be 1.5% and 2% [25], respectively.
Systematic uncertainty from ∆E requirement is studied using the control samples of D 0 → K − π + π 0 and D + → K + π + π − for the signal processes with and without π 0 in final states, where the control samples are se- lected with the ST method. We set [µ − 3.5σ, µ + 2.5σ] as a nominal ∆E window for the signal, where µ and σ are the mean and width values of ∆E obtained by fitting. Then we vary ∆E window by 0.5σ for both sides, and the resulting difference on the change of efficiency between data and MC simulation is taken as the systematic uncertainty.
To study the systematic uncertainty associated with FSR recovery process, we also obtain the alternative detection efficiency without the FSR recovery process, and the change on the efficiency is taken as the systematic uncertainty.
The difference on the detector efficiency between the one obtained with the phase space generator, and the average one with the m N -dependent case, is taken as the systematic uncertainty associated with the modeling for detector efficiency.
Systematic uncertainty associated with the fitting on the M BC distribution arises from fitting range, signal shape, and background shape. We performed the alternative fits by varying the fitting range from [1.84, 1.89] to [1.845, 1.89] GeV/c 2 , the width of convolved Gaussian for signal shape within one standard deviation, and the background shape from ARGUS function to the cocktail MC simulated shape, individually. The relative changes on the signal yields will be taken as the corresponding systematic uncertainties, and are found to be negligible comparing to statistical fluctuations.
All the discussed systematic uncertainties are summarized in Table II. Assuming they are independent, the total systematic uncertainty is the quadrature sum of individual ones.

VI. RESULT AND DISCUSSION
A. Upper limits for D → Kπe + e + decays Taking into account the effect of systematic uncertainties, we calculate the ULs on the BFs for the LNV ∆L = 2 decays D 0 → K − π − e + e + , D + → K 0 S π − e + e + and D + → K − π 0 e + e + individually according to Eq. (2) based on the Bayesian method [26]. A series fits on the M BC distribution are carried out by fixing the BF at different values, and the resultant curve of likelihood values L as a function of the BF is convolved with a Gaussian function, which is of the width of the overall systematic uncertainty, and is normalized to the maximum value of 1. The ULs on the BF at the 90% CL, B UL sig , for the different processes are the values that yield 90% of the likelihood integral over BF from zero to infinity, as listed in Table III.

B. Searching for Majorana neutrino
With the above three decay processes, the Majorana neutrino can be searched for by studying the decay chains D 0 → K − e + ν N (π − e + ), D + → K 0 S e + ν N (π − e + ) or D + → π 0 e + ν N (K − e + ), a narrow peak will be present on the distribution of π − e + (K − e + ) invariant mass if signal exists. Comparing to the other two decay channels, the D + → π 0 e + ν N (K − e + ) is expected to be suppressed by a factor of 1/20 because of the smaller CKM factors. So in this analysis, the Majorana neutrino is searched for in the processes D 0 → K − e + ν N (π − e + ) and D + → K 0 S e + ν N (π − e + ) with different mass hypotheses, i.e., between 0.25 to 1.0 GeV/c 2 with an interval of 0.05 GeV/c 2 .
Based on the above selection criteria, to search for the Majorana neutrino with a given mass, m N , the candidate events are selected by further requiring the invariant mass of any π − e + combination (two e + per event), M π − e + , within the range of [m N − 3σ, m N + 3σ], where σ is the resolution of the M π − e + distribution obtained by studying the signal MC sample. The number counting method is used to extract the signal yields due to very little events survived. We count the number of signal candidates within the M BC region of [1.859, 1.872] ([1.865, 1.875]) GeV/c 2 for the decay D 0 → K − e + ν N (π − e + ) TABLE II. Summary of the relative systematic uncertainties for D → Kπe + e + processes (in percent). Here '-' denotes that there is no corresponding systematic uncertainty, and '...' means that the corresponding systematic uncertainty is negligible.

Source
Relative systematic uncertainty (%)  Fig. 2. The ULs on the BFs of Majorana neutrino case are calculated through the profile likelihood method incorporating the systematic uncertainty by implementing the package trolke [27,28] in root framework, where the numbers of events in the signal and sideband regions are assumed to be the Poisson distributions, and the efficiency is Gaussian distribution. The ULs on the BFs at the 90% CL as a function of m N is at the level of 10 −7 ∼ 10 −6 , as shown in Fig. 3.
Based on the measured BFs, the mixing matrix element |V eN | 2 of a positron with the heavy Majorana neutrino in the charged current interaction [9,10] as a function of m N can be obtained by Eq. (3) [12], where the decay width Γ(m N , V eN (m N )) is proportional to its BF, and Γ(m N , V ′ eN (m N )) is related to the BF given in Tables 4 and 5 of Ref. [12], based on the assumptions that the Majorana neutrino is on-shell and its width is negligible compared to the neutrino mass. The mixing matrix element |V ′ eN (m N )| 2 is derived from a reanalysis of neutrinoless double beta decay experimental data [29]. The resultant ULs on the mixing matrix element |V eN | 2 as a function of m N , which are also depicted in Fig. 3, provide the additional and complementary information about the bounds on the |V eN | 2 in D meson decays.

VII. SUMMARY
Using the data sample with the integral luminosity of 2.93 fb −1 collected at the C.M. energy √ s = 3.773 GeV, we perform a search for LNV ∆L = 2 decays of D → Kπe + e + as well as search for Majorana neutrino with different mass hypotheses. No evidence of signal is found. Therefore, by applying the Bayesian approach, we place the 90% CL the ULs on the decay BF for the D 0 → K − π − e + e + , D + → K 0 S π − e + e + and D + → K − π 0 e + e + to be 2.7 × 10 −6 , 3.3 × 10 −6 and 8.5×10 −6 , respectively. We also measured the ULs on the BF at the 90% CL for the decays D 0 → K − e + ν N (π − e + ) and D + → K 0 S e + ν N (π − e + ) with different m N hypotheses within the range between 0.25 and 1.0 GeV/c 2 , which UL on BF at the 90% CL are of the level 10 −7 ∼ 10 −6 . The constraints on the mixing element |V eN | 2 depending on m N is also evaluated based on the measured BFs. The presenting results provide the supplementary information to study the mixing between the heavy Majorana neutrino and the standard model neutrino ν e of flavour e in D meson decays.