Search for rare or forbidden decays of the $D^{0}$ meson

We present a search for nine lepton-number-violating and three lepton-flavor-violating neutral charm decays of the type $D^0\rightarrow h^{\prime -} h^{-}\ell^{\prime +} \ell^{+}$ and $D^0\rightarrow h^{\prime -} h^{+}\ell^{\prime\pm} \ell^{\mp}$, where $h$ and $h^{\prime}$ represent a $K$ or $\pi$ meson and $\ell$ and $\ell^{\prime}$ an electron or muon. The analysis is based on $468$ fb$^{-1}$ of $e^+e^-$ annihilation data collected at or close to the $Y(4S)$ resonance with the BaBar detector at the SLAC National Accelerator Laboratory. No significant signal is observed for any of the twelve modes and we establish 90% confidence level upper limits on the branching fractions in the range $(1.0 - 30.6)\times 10^{-7}$. The limits are between one and three orders of magnitude times more stringent than previous measurements.

Accelerator Laboratory. No significant signal is observed for any of the twelve modes, and we establish 90% confidence level upper limits on the branching fractions in the range (1.0−30.6)×10 −7 . The limits are between 1 and 3 orders of magnitude more stringent than previous measurements.
PACS numbers: 13.25.Ft,11.30.Fs Lepton-flavor-violating and lepton-number-violating neutral charm decays can be used to investigate physics beyond the standard model (SM) of particle physics. A potential set of decays for study are of the form D 0 → h − h − + + and D 0 → h − h + ± ∓ , where h and h represent a K or π meson and and an electron or muon [1].
The D 0 → h − h + ± ∓ decay modes with two opposite-charge, different-flavor leptons in the final state are lepton-flavor-violating (LFV). They are essentially prohibited in the SM because they can occur only through lepton mixing [2]. The D 0 → h − h − + + decay modes with two same-charge leptons are both leptonflavor violating and lepton-number violating (LNV) and are forbidden in the SM in low-energy collisions or decays. However, LNV processes can occur in extremely high-energy or high-density interactions [3].
Lepton-number violation is a necessary condition for leptogenesis as an explanation of the baryon asymmetry of the Universe [4]. If neutrinos are Majorana fermions, the neutrino and antineutrino are the same particle, and some LNV processes become possible [5]. Many models beyond the SM allow lepton-number violation. Most models have made predictions for, or used constraints from, three-body decays of the form D → M l l or B → M l l, where M is a meson [6][7][8][9][10][11][12]. However, some models that consider LFV and LNV four-body charm decays predict branching fractions up to O(10 −6 ) to O(10 −5 ), approaching those accessible with current data [11][12][13].
In this report we present a search for nine D 0 → h − h − + + LNV decays and three D 0 → h − h + ± ∓ LFV decays, with data recorded with the BABAR detector at the PEP-II asymmetric-energy e + e − collider operated at the SLAC National Accelerator Laboratory. The data sample corresponds to 424 fb −1 of e + e − collisions collected at the center-of-mass (c.m.) energy of the Υ (4S) resonance (on peak) and an additional 44 fb −1 of data collected 40 MeV below the Υ (4S) resonance (off peak) [24]. The branching fractions for signal modes with zero, one, or two kaons in the final state are measured relative to the normalization decays D 0 → π − π + π + π − , D 0 → K − π + π + π − , and D 0 → K − K + π + π − , respectively. The D 0 mesons are identified from the decay D * + → D 0 π + produced in e + e − → cc events.
The BABAR detector is described in detail in Refs. [25,26]. Charged particles are reconstructed as tracks with a five-layer silicon vertex detector and a 40-layer drift chamber inside a 1.5 T solenoidal magnet. An electromagnetic calorimeter comprised of 6580 CsI(Tl) crystals is used to identify electrons and photons. A ringimaging Cherenkov detector is used to identify charged hadrons and to provide additional lepton identification information. Muons are identified with an instrumented magnetic-flux return.
Monte Carlo (MC) simulation is used to investigate sources of background contamination, evaluate selection efficiencies, cross-check the selection procedure, and for studies of systematic effects. The signal and normalization channels are simulated with the EvtGen package [27]. We generate the signal channel decays uniformly throughout the four-body phase space, while the normalization modes include two-body and three-body intermediate resonances, as well as nonresonant decays. We also generate e + e − → qq (q = u, d, s, c), dimuon, Bhabha elastic e + e − scattering, BB background, and two-photon events [28,29]. Final-state radiation is generated using Photos [30]. The detector response is simulated with GEANT 4 [31,32].
In order to optimize the event reconstruction, candidate selection criteria, multivariate analysis training, and fit procedure, a rectangular area in the m(D 0 ) versus ∆m = m(D * + ) − m(D 0 ) plane is defined, where m(D * + ) and m(D 0 ) are the reconstructed masses of the D * + and D 0 candidates, respectively. This optimization region is kept hidden (blinded) in data until the analysis steps are finalized. The blinded region is approximately 3 times the width of the ∆m and m(D 0 ) resolutions. The ∆m region is 0.1447 < ∆m < 0.1462 GeV/c 2 for all modes. The m(D 0 ) signal peak distribution is asymmetric due to bremsstrahlung emission. The upper m(D 0 ) bound on the blinded region is 1.874 GeV/c 2 for all modes, and the lower bound is 1.848 GeV/c 2 , 1.852 GeV/c 2 , and 1.856 GeV/c 2 for modes with two, one or no electrons, respectively.
Events are required to contain at least five charged tracks. Particle identification (PID) criteria are applied to all the charged tracks to identify kaons, pions, electrons, and muons [26,33]. For modes with two kaons in the final state, the PID requirement on the kaons is relaxed compared to the single-kaon modes. This increases the reconstruction efficiency for the modes with two kaons, with little increase in backgrounds or misidentified candidates. The PID efficiency depends on the track momentum, and is in the range 0.96 − 0.99 for electrons, 0.60 − 0.95 for muons, and 0.90 − 0.98 for kaons and pions. The misidentification probability is typically less than 0.03 for all selection criteria, except for the pion selection criteria, where the muon misidentification rate can be as high as 0.35 at low momentum.
Candidate D 0 mesons are formed from four charged tracks reconstructed with the appropriate mass hypotheses for the signal and normalization decays. The four tracks must form a good-quality vertex with a χ 2 probability for the vertex fit greater than 0.005. A bremsstrahlung energy recovery algorithm is applied to electrons [16]. The invariant mass of any e + e − pair is required to be greater than 0.1 GeV/c 2 . For the normalization modes, the reconstructed D 0 meson mass is required to be in the range 1.81 < m(D 0 ) < 1.91 GeV/c 2 , while for the signal modes, m(D 0 ) must be in the blinded m(D 0 ) range defined above.
The candidate D * + is formed by combining the D 0 candidate with a charged pion with a momentum in the laboratory frame greater than 0.1 GeV/c. For the normalization mode D 0 → K − π + π + π − , this pion is required to have a charge opposite to that of the finalstate kaon. A vertex fit is performed with the D 0 mass constrained to its known value [20] and the requirement that the D 0 meson and the pion originate from the PEP-II interaction region. The χ 2 probability of the fit is required to be greater than 0.005. For signal modes with two kaons, the mass difference ∆m is required to be 0.141 < ∆m < 0.201 GeV/c 2 . Signal modes with fewer than two kaons have almost no candidates beyond ∆m = 0.149 GeV/c 2 , and the range for these modes is restricted to 0.141 < ∆m < 0.149 GeV/c 2 .
After the application of the D * + vertex fit, the D 0 candidate momentum in the c.m. system, p * , is required to be greater than 2.4 GeV/c. This removes most sources of combinatorial background and also charm hadrons produced in B decays, which are limited to p * < ∼ 2.2 GeV/c [34]. Remaining backgrounds are mainly radiative Bhabha scattering, initial-state radiation, and two-photon events, which are all rich in electrons and positrons. We suppress these backgrounds by requiring that the PID signatures of the hadron candidates be inconsistent with the electron hypothesis.
Hadronic D 0 decays with large branching fractions, where one or more charged tracks are misidentified as leptons, will usually have reconstructed D 0 masses well away from the known D 0 mass [20]. To ensure rejection of this type of background. the D 0 candidate is also reconstructed assuming the kaon or pion mass hypothesis for the lepton candidates. If the resulting D 0 candidate mass is within 20 MeV/c 2 of the known D 0 mass, and if |∆m| < 2 MeV/c 2 , the event is discarded. After these criteria are applied, the background from these hadronic decays is negligible.
Two particular sources of background are semileptonic charm decays in which a charged hadron is misidentified as a lepton; and charm decays in which the final state contains a neutral particle or more than four charged tracks. In both cases, tracks can be selected from elsewhere in the event to form a D 0 candidate. To reject these backgrounds, a multivariate selection based on a Fisher discriminant is applied [35]. The discriminant uses nine input observables: the momenta of the four tracks used to form the D 0 candidate; the thrust and sphericity of the D * + candidate [36]; the angle between the D * + meson candidate sphericity axis and the sphericity axis defined by the charged particles in the rest of the event (ROE); the angle between the D * + meson candidate thrust axis and the thrust axis defined by the charged particles in the ROE; and the second Fox-Wolfram moment [37] calculated from the entire event using both charged and neutral particles. The input observables are determined in the laboratory frame after the application of the D * + vertex fit. The discriminant is trained and tested using MC for the signal modes; for the background, data outside the optimization region, together with e + e − → cc MC samples, are used. The training is performed independently for each signal mode. A requirement on the Fisher discriminant output is chosen such that approximately 90% of the simulation signal candidates are accepted. Depending on the signal mode, this rejects 30% to 50% of the background in data.
The cross feed to one signal mode from the other eleven is estimated from MC samples to be < ∼ 0.5% in all cases, assuming equal branching fractions for all signal modes. The cross feed to a specific normalization mode from the other two normalization modes is predicted from simulation to be < ∼ 0.7%, where the branching fractions are taken from Ref. [20]. Multiple candidates occur in 4.5% to 7.1% of simulated signal events and in 2.4% to 4.4% of the normalization events in data. If two or more candidates are found in an event, the one with the highest D * + vertex χ 2 probability is selected. After the application of all selection criteria and corrections for small differences between data and MC simulation in tracking and PID performance derived from high purity control samples [26], the reconstruction efficiency sig for the simulated signal decays is between 3.2% and 6.2%, depending on the mode. For the normalization decays, the reconstruction efficiency norm is between 19.2% and 24.7%. The difference between sig and norm is mainly due to the momentum dependence of the lepton PID [26].
The signal mode branching fraction B sig is determined relative to that of the normalization decay using where B norm is the normalization mode branching fraction [20], and N sig and N norm are the fitted yields of the signal and normalization mode decays, respectively. The symbols L sig and L norm represent the integrated luminosities of the data samples used for the signal (468.2 ± 2.0 fb −1 ) and the normalization decays (39.3 ± 0.2 fb −1 ), respectively [24]. For the signal modes, we use the on-peak and off-peak data samples, while the normalization modes use only a subset of the off-peak data. Each normalization mode yield N norm is extracted by performing an extended two-dimensional unbinned maximum likelihood (ML) fit [38] to the observables ∆m and m(D 0 ) in the range 0.141 < ∆m < 0.149 GeV/c 2 and 1.81 < m(D 0 ) < 1.91 GeV/c 2 . The measured ∆m and m(D 0 ) values are not correlated and are treated as independent observables in the fits. The probability density functions (PDFs) in the fits depend on the normalization mode and use sums of multiple Cruijff [16] and Crystal Ball [39] functions in both ∆m and m(D 0 ). The functions for each observable use a common mean. The background is modeled with an ARGUS threshold function [40] for ∆m and a Chebyshev polynomial for m(D 0 ). The ARGUS end point parameter is fixed to the kinematic threshold for a D * + → D 0 π + decay. All other PDF parameters, together with the normalization mode and background yields, are allowed to vary in the fit. The fitted yields and reconstruction efficiencies for the normalization modes are given in Table I.   TABLE I. Summary of fitted candidate yields with statistical uncertainties, systematic uncertainties, and reconstruction efficiencies for the three normalization modes.

Decay mode
Nnorm Systematic norm (candidates) (%) (%) D 0 → K − π + π + π − 260 870 ± 520 4.7 20.1 ± 0.2 D 0 → K − K + π + π − 8480 ± 110 6.6 19.2 ± 0.2 D 0 → π − π + π + π − 28 470 ± 220 6.8 24.7 ± 0.2 Each signal mode yield N sig is extracted by performing the ML fit with the single observable ∆m in the range 0.141 < ∆m < 0.201 GeV/c 2 for signal modes with two kaons and 0.141 < ∆m < 0.149 GeV/c 2 for all other signal modes. The signal PDF is a Cruijff function with parameters obtained by fitting the signal MC. The background is modeled with an ARGUS function with an end point that is set to the same value that is used for the normalization modes. The signal PDF parameters and the end point parameter are fixed in the fit. All other background parameters and the signal and background yields are allowed to vary. Figures 1 and 2 show the results of the fits to the ∆m distributions for the twelve signal modes.  We test the performance of the ML fit for the normalization modes by generating ensembles of MC samples from the normalization and background PDF distributions. The mean numbers of normalization and background candidates used in the ensembles are taken from the fits to the data. The numbers of background and normalization mode candidates are allowed to fluctuate according to a Poisson distribution and all background and normalization mode PDF parameters are allowed to vary. No significant biases are observed in fitted yields of the normalization modes. The same procedure is repeated for the ML fit to signal modes, with ensembles of MC samples generated from the background PDF distributions only, assuming a signal yield of zero. The signal PDF parameters are fixed to the values used for the fits to the data but the signal yield is allowed to vary. The biases in the fitted signal yields are less than ±0.2 for all modes, and these are subtracted from the fitted yields before calculating the signal branching fractions.
The main sources of systematic uncertainties in the signal yields are associated with the model parametrizations used in the fits to the signal modes and backgrounds, the fit biases, and the limited MC and data sample sizes available for the optimization of the Fisher discriminants.
The uncertainties associated with the fit model parametrizations of the signal modes are estimated by repeating the fits with alternative PDFs. This involves swapping the Cruijff and Crystal Ball functions, using Gaussian functions with different asymmetric widths, and changing the number of functions used. For the background, the order of the polynomials is changed and the ARGUS function is replaced by a second-order polynomial. The fits are also performed with the fixed signal parameters allowed to vary within the statistical uncertainties obtained from fits to the signal MC samples. The systematic uncertainty is taken as half the maximum deviation from the default fit.
The systematic uncertainties in the corrections on the fit biases for the signal yields are taken to be the statistical uncertainties on the ensembles of fits to the MC samples described above. The systematic uncertainty due to knowledge of the Fisher discriminant shape is obtained by varying the value of the selection criterion for the Fisher discriminant, changing the size of the blinded optimization region, and retraining the Fisher discriminant using a training sample with a different set of MC samples. The uncertainty is taken as half the maximum difference from the yield obtained with the default Fisher discriminant criterion. Summed together, the total systematic uncertainties in the signal yield are between 0.4 and 1.9 events, depending on the mode.
Systematic uncertainties that impact the calculation of the branching fractions of the signal modes are due to assumptions made about the distributions of the finalstate particles in the signal simulation modeling, the model parametrizations used in the fits to the normalization modes, the normalization mode branching fractions, tracking and PID efficiencies, and luminosity.
Since the decay mechanism of the signal modes is unknown, we vary the angular distributions of the simulated final-state particles from the D 0 signal decay, where the three angular variables are defined following the prescription of Ref. [41]. We weight the reconstruction efficiencies of the phase-space simulation samples as a function of the angular-variable distributions, trying combinations of sin, cos, sin 2 , and cos 2 functions. Half the maximum change in the average reconstruction efficiency is assigned as a systematic uncertainty.
Uncertainties associated with the fit model parametrizations of the normalization modes are estimated by repeating the fits with alternative PDFs for the normalization modes and backgrounds. Uncertainties in the normalization mode branching fractions are taken from Ref. [20]. We include reconstruction efficiency uncertainties of 0.8% per track for the leptons and 0.7% for the kaon and pion [42]. For the PID efficiencies, we assign an uncertainty of 0.7% per track for electrons, 1.0% for muons, 0.2% for pions, and 1.1% for kaons [26]. A systematic uncertainty of 0.43% is associated with our knowledge of the luminosities L norm and L sig [24]. The total systematic uncertainties in the signal efficiencies are between 5% and 19%, depending on the mode.
We use the frequentist approach of Feldman and Cousins [43] to determine 90% C.L. bands that relate the true values of the branching fractions to the measured signal yields. When computing the limits, the systematic uncertainties are combined in quadrature with the statistical uncertainties in the fitted signal yields.
The signal yields for all the signal modes are compatible with zero. Table II gives the fitted signal yields, reconstruction efficiencies, branching fractions with statistical and systematic uncertainties, and 90% C.L. branching fraction upper limits for the signal modes.
In summary, we report 90% C.L. upper limits on the branching fractions for nine lepton-numberviolating D 0 → h − h − + + decays and three leptonflavor-violating D 0 → h − h + ± ∓ decays. The analysis TABLE II. Summary of fitted signal yields with statistical and systematic uncertainties, reconstruction efficiencies, branching fractions with statistical and systematic uncertainties, 90% C.L. branching fraction upper limits (U.L.), and the previous limits [20,21]. The branching fraction systematic uncertainties take into account correlations and cancellations between the signal and normalization modes and include the uncertainties in the normalization mode branching fractions. is based on a sample of e + e − annihilation data collected with the BABAR detector, corresponding to an integrated luminosity of 468.2±2.0 fb −1 . The limits are in the range (1.0 − 30.6) × 10 −7 and are between 1 and 3 orders of magnitude more stringent than previous results.