Observation of the decay $D^0\rightarrow K^-\pi^+e^+e^-$

We report the observation of the rare charm decay $D^0\rightarrow K^-\pi^+e^+e^-$, based on $468$ fb$^{-1}$ of $e^+e^-$ annihilation data collected at or close to the center-of-mass energy of the $Y(4S)$ resonance with the BaBar detector at the SLAC National Accelerator Laboratory. We find the branching fraction in the invariant mass range $0.675<m(e^+e^-)<0.875 GeV/c^2$ of the electron-positron pair to be ${\cal B}(D^0\rightarrow K^-\pi^+e^+e^-) = (4.0\pm0.5\pm0.2\pm0.1) \times 10^{-6}$, where the first uncertainty is statistical, the second systematic, and the third due to the uncertainty in the branching fraction of the decay $D^0\rightarrow K^-\pi^+\pi^+\pi^-$ used as a normalization mode. The significance of the observation corresponds to 9.7 standard deviations including systematic uncertainties. This result is consistent with the recently reported $D^0\rightarrow K^-\pi^+\mu^+\mu^-$ branching fraction, measured in the same invariant mass range, and with the value expected in the Standard Model. In a set of regions of $m(e^+e^-)$ where long-distance effects are potentially small, we determine a 90% confidence level upper limit on the branching fraction ${\cal B}(D^0\rightarrow K^-\pi^+e^+e^-)<3.1 \times 10^{-6}$.

systematic, and the third due to the uncertainty in the branching fraction of the decay D 0 → K − π + π + π − used as a normalization mode. The significance of the observation corresponds to 9.7 standard deviations including systematic uncertainties. This result is consistent with the recently reported D 0 → K − π + µ + µ − branching fraction, measured in the same invariant mass range.
PACS numbers: 13.20.Fc,11.30.Hv The decay D 0 → K − π + e + e − [1] is an example of a flavor-changing neutral-current (FCNC) process, which is expected to be very rare in the standard model (SM) as it cannot occur at tree level and is suppressed by the Glashow-Iliopoulos-Maiani (GIM) mechanism at loop level [2]. Short-distance contributions to the D 0 → K − π + e + e − branching fraction proceed through loop and box diagrams [3] and are expected to be O(10 −9 ). However, decays with long-distance contributions, such D 0 → XV involving a vector meson, V, decaying to two leptons, could contribute at the level of O(10 −6 ) [3][4][5][6][7].
Certain physics models beyond the standard model, such as minimal supersymmetric or R-parity-violating supersymmetric theories, predict branching fractions as high as O(10 −5 ) [3,[7][8][9][10]. As virtual particles can enter in the one-loop processes, this type of decay can be used to study new physics processes at large mass scales. These processes could potentially be detected in regions where intermediate vector mesons do not dominate.
Over the last few years there have been a number of measurements of the decays of B mesons to final states involving one or more charged leptons. Some of these measurements suggest a possible deviation from the assumption that all leptons couple equally [11][12][13][14][15][16][17][18][19][20]. The possibility therefore exists that a deviation from lepton universality will be seen in D meson decays.
We report here the observation of the decay D 0 → K − π + e + e − 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 energy of the Υ (4S) resonance (onpeak) and an additional 44 fb −1 of data collected 40 MeV below the Υ (4S) resonance (offpeak) [25]. The signal branching fraction B(D 0 → K − π + e + e − ) is measured relative to the normalization decay D 0 → K − π + π + π − . The D 0 mesons are reconstructed from the decay D * + → D 0 π + produced in e + e − → cc events. The use of this decay chain increases the purity of the sample at the expense of a smaller number of reconstructed D 0 mesons.
The BABAR detector is described in detail in Refs. [26,27]. 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 evaluate the level of background contamination and selection efficiencies. Simulated events are also used to cross-check the selection procedure and for studies of systematic effects. The signal and normalization channels are simulated with the EvtGen package [28]. We generate the signal channel decay with a phase-space model, while the normalization mode includes two-body and threebody intermediate resonances, as well as nonresonant decays. For background studies, we generate e + e − → qq (q = u, d, s, c), dimuon, Bhabha elastic e + e − scattering, BB background, and two-photon events [29,30]. The background samples are produced with an integrated luminosity approximately six times that of the data. Finalstate radiation is provided by PHOTOS [31]. The detector response is simulated with GEANT4 [32,33]. All simulated events are reconstructed in the same manner as the data.
Events are required to contain at least five charged tracks. Candidate D 0 mesons are formed from four charged tracks reconstructed with the appropriate mass hypothesis for the D 0 → K − π + e + e − and D 0 → K − π + π + π − decays. Particle identification (PID) is applied to the charged tracks and the same criteria are applied to the signal and normalization modes [27,34]. The four tracks must form a good-quality vertex with a χ 2 probability for the vertex fit greater than 0.005. In the case of D 0 → K − π + e + e − , a bremsstrahlung energy recovery algorithm is applied to the electrons, in which the energy of photon showers that are within a small angle (typically 35 mrad) of the initial electron direction are added to the energy of the electron candidate. The electron-positron pair must have an invariant mass m(e + e − ) > 0.1 GeV/c 2 . The D 0 candidate momentum in the PEP-II center-of-mass system, p * , must be greater than 2.4 GeV/c. The requirement for five charged tracks strongly suppresses backgrounds from QED processes. The p * criterion removes most sources of combinatorial background and also charm hadrons produced in B decays, which are kinematically limited to less than ∼2.2 GeV/c.
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. The pion is required to have a charge opposite to that of the kaon in the D 0 decay. A vertex fit is performed with the D 0 mass constrained to its known value and the requirement that the D 0 meson and the pion originate from the interaction region. The χ 2 probability of the fit is required to be greater than 0.005. The D 0 meson mass m(D 0 ) must be in the range 1.81 < m(D 0 ) < 1.91 GeV/c 2 and the mass difference, ∆m = m(D * + ) − m(D 0 ), between the reconstructed masses of the D * + and D 0 candidates is required to satisfy 0.143 < ∆m < 0.148 GeV/c 2 .
To reject misreconstructed D 0 → K − π + e + e − candidates that originate from D 0 hadronic decays with large branching fractions where one or more charged tracks are misidentified by the PID, the candidate is reconstructed assuming the kaon or pion mass hypothesis for the leptons. If the resulting candidate m(D 0 ) is within 20 MeV/c 2 of the known D 0 mass [35] and |∆m| < 2 MeV/c 2 , the event is discarded. After these criteria are applied, the background from these hadronic decays is negligible. Multiple candidates occur in less than 4% of simulated D 0 → K − π + π + π − decays and in less than 2% of D 0 → K − π + e + e − decays. If two or more candidates are found in an event, the one with the highest vertex χ 2 probability is selected.
After the application of all selection criteria and corrections for small differences between data and MC tracking and PID performance, the average reconstruction efficiency for the D 0 → K − π + π + π − decay isˆ norm = (20.1±0.2)%. For the D 0 → K − π + e + e − decay, the average reconstruction efficiencyˆ sig varies between 5.0% and 8.9% depending on the m(e + e − ) mass range, as shown in Table I. The remaining background comes predominantly from e + e − → cc events. No evidence is found in MC for backgrounds that peak in the m(D 0 ) and ∆m mass range.
The D 0 → K − π + e + e − branching fraction is determined relative to that of the normalization decay channel where B(D 0 → K − π + π + π − ) is the branching fraction of the normalization mode [35], and N norm andˆ norm are the D 0 → K − π + π + π − fitted yield and the reconstruction efficiency calculated from simulated D 0 → K − π + π + π − decays, respectively. The fitted D 0 → K − π + e + e − signal yield is represented by N sig and i sig is the reconstruction efficiency for each signal candidate i calculated as a function of m(e + e − ) and m(K − π + ) from simulated D 0 → K − π + e + e − decays. The symbols L sig and L norm represent the integrated luminosities used for the signal D 0 → K − π + e + e − decay (468.2 ± 2.0 fb −1 ) and the normalization D 0 → K − π + π + π − decay (39.3 ± 0.2 fb −1 ), respectively [25]. The signal mode uses both the onpeak and offpeak data samples while the normalization mode uses only a subset of the offpeak data.
The D 0 → K − π + e + e − and D 0 → K − π + π + π − yields are determined from extended unbinned maximum likelihood fits to the ∆m and the four-body mass distributions. The ∆m and the four-body mass distributions are not correlated and are treated as independent observables in the fit. For the D 0 → K − π + e + e − signal, a Gaussianlike function with different lower and upper widths is used for both ∆m and m(K − π + e + e − ). This asymmetric function is used in order to describe bremsstrahlung emission of photons by the electrons. The background in the D 0 → K − π + e + e − channel is modeled with an ARGUS threshold function [36] for ∆m and a firstorder Chebychev polynomial for m(K − π + e + e − ). For the D 0 → K − π + π + π − normalization mode, the ∆m and m(K − π + π − π + ) distributions are each represented by two Cruijff functions with shared means [37]. The background is represented by an ARGUS threshold function for ∆m and a second-order Chebychev polynomial for m(K − π + π − π + ). All yields and shape parameters are allowed to vary in the fits except for the ARGUS function threshold end point, which is set to the kinematic threshold for a D * + → D 0 π + decay.
The fitted yield for the D 0 → K − π + π + π − normalization data sample is 260 870 ± 520 candidates. For the D 0 → K − π + e + e − signal mode, the fitted yield is 68 ± 9 candidates in the range 0.675 < m(e + e − ) < 0.875 GeV/c 2 . The significance S = √ −2∆ ln L of the signal yield in this mass range, including statistical and systematic uncertainties, is 9.7 standard deviations (σ), where ∆ ln L is the change in the log-likelihood from the maximum value to the value when the number of D 0 → K − π + e + e − signal candidates is set to N sig = 0. Figure 1 shows the results of the fit to the m(K − π + e + e − ) and ∆m distributions of the D 0 → K − π + e + e − signal mode in the mass range 0.675 < m(e + e − ) < 0.875 GeV/c 2 . Figure 2 shows the projection of the fit to the D 0 → K − π + e + e − signal mode as a function of m(e + e − ) and m(K − π + ), where the background has been subtracted using the sPlot technique [38]. A peaking structure is visible in m(e + e − ) centered near the ρ 0 mass. A broader structure is seen in m(K − π + ) near the known mass of the K * (892) 0 meson. Both distributions are similar to the distributions shown in Ref. [21] for the decay D 0 → K − π + µ + µ − . Figure 3 shows the m(e + e − ) and m(K − π + ) distributions for the expanded mass region m(e + e − ) > 0.1 GeV/c 2 , where the back-ground has been subtracted using the sPlot technique. A peak can be seen near the mass of the π 0 , compatible with the decay D 0 → K − π + π 0 , π 0 → e + e − γ, where the γ has not been reconstructed. There are indications in Fig. 3(b) of other intermediate states, such as φ and η mesons, but they are not statistically significant.  We test the performance of the maximum likelihood fit by generating ensembles of MC samples from both the PDF distributions and the fully simulated MC events. The mean number of signal, normalization, and background yields used in the ensembles is taken from the fits to the data sample. The yields are allowed to fluctuate according to a Poisson distribution and all fit parameters are allowed to vary. No significant bias is observed in the normalization mode. The largest fit bias observed in the signal mode is 0.4 ± 0.1 candidates. The biases are much smaller than the statistical uncertainties in the yields. The fit biases are subtracted from the fitted yields before calculating the signal branching fractions.
To cross-check the normalization procedure, the signal mode D 0 → K − π + e + e − in Eq. (1) is replaced with the decay D 0 → K − π + , which has a well-known branching fraction [35]. We use the onpeak data sample to reconstruct the D 0 → K − π + channel. The D 0 → K − π + decay is selected using the same criteria as used for the D 0 → K − π + π + π − mode. The D 0 → K − π + yield is determined using an unbinned maximum likelihood fit to ∆m and the two-body invariant mass m(K − π + ). Three Crystal Ball functions [39] with shared means are used for the D 0 → K − π + signal ∆m and m(K − π + ) distributions. The backgrounds are represented by an ARGUS function for ∆m and a second-order Chebychev polynomial for m(K − π + ). The D 0 → K − π + fitted signal yield is 1 881 950 ± 1380 candidates with an average reconstruction efficiency ofˆ sig = (27.4 ± 0.2)%. We determine B(D 0 → K − π + ) = (3.98 ± 0.08 ± 0.10)%, where the uncertainties are statistical and systematic, respectively; the current world-average is (3.89 ± 0.04)% [35]. Similar compatibility with the B(D 0 → K − π + ) worldaverage, but with larger uncertainties, is achieved when using the four-body decay modes D 0 → K + K − π + π − and D 0 → π − π + π + π − instead of D 0 → K − π + π + π − for the normalization mode in Eq. (1). The main sources of systematic uncertainty are associated with the model parameterizations used in the fits and the normalization procedure, signal MC model, fit bias, tracking and PID efficiencies, luminosity, and the normalization mode branching fraction. Some of the tracking and PID systematic effects cancel in the branching fraction determination since they affect both the signal and normalization modes.
Systematic uncertainties associated with the model parameterization are estimated by repeating the fit with the D 0 → K − π + e + e − signal parameters for the ∆m and four-body distributions fixed to values taken from simulation. Alternative fits are also performed with the default functions for the signal and normalization modes replaced with alternative functions. The resulting uncertainties are 1.9% and 1.0% for the signal and normalization yields, respectively.
In the mass range 0.675 < m(e + e − ) < 0.875 GeV/c 2 , we replace the signal phase-space simulation model with a model assuming D 0 → K * (892) 0 ρ 0 with K * (892) 0 → K − π + and ρ 0 → e + e − and assign half the difference with the default reconstruction efficiency as a systematic uncertainty, equivalent to a relative change of 1.8%. We also use this number as an estimate of the relative change in other regions of m(e + e − ) and m(K − π + ) where no suitable alternative simulation model exists.
The systematic uncertainty in the fit bias for the signal yield is taken from the ensemble of fits to the MC data samples and we attribute a value of half the largest fit bias found, ±0.2 candidates. To account for imperfect knowledge of the tracking efficiency, we assign an uncertainty of 0.8% per track for the leptons and 0.7% for the kaon and pion [40]. For the PID, we estimate an uncertainty of 0.7% per electron, 0.2% per pion, and 1.1% per kaon [27]. A systematic uncertainty of 0.8% is associated with the knowledge of the luminosity ratio, L norm /L sig [25].
The overall systematic uncertainty in the yields is 5.3% for the signal and 3.6% for the normalization mode. As the PID and tracking systematic uncertainties of the kaons and pions are correlated and cancel, the combined systematic uncertainty in the D 0 → K − π + e + e − branching fraction is 3.8%, where the uncertainty in the D 0 → K − π + π + π − branching fraction is excluded [35]. Table I summarizes the fitted signal yields and branching fractions for different m(e + e − ) mass regions. The branching fraction B(D 0 → K − π + e + e − ) in the mass range 0.675 < m(e + e − ) < 0.875 GeV/c 2 is determined to be (4.0±0.5±0.2±0.1)×10 −6 , where the first uncertainty is statistical, the second systematic, and the third comes from the uncertainty in B(D 0 → K − π + π + π − ) [35]. This result is compatible within the uncertainties with B(D 0 → K − π + µ + µ − ) reported in Ref. [21]. In the mass range m(e + e − ) > 0.2 GeV/c 2 , but excluding 0.675 < m(e + e − ) < 0.875 GeV/c 2 , B(D 0 → K − π + e + e − ) = (4.8 ± 0.7 ± 0.3 ± 0.1) × 10 −6 . Although there is no statistically significant evidence for individual intermediate vector mesons in this m(e + e − ) range, long-distance effects are not ruled out. We do not report branching fractions for m(e + e − ) mass regions that include the potential intermediate decay π 0 → e + e − γ as the photon is not explicitly reconstructed.
In summary, we have presented the first observation of the decay D 0 → K − π + e + e − . The analysis is based on a sample of e + e − annihilation data collected with the BABAR detector, corresponding to an integrated luminosity of 468 fb −1 . The branching fraction in the mass range 0.675 < m(e + e − ) < 0.875 GeV/c 2 is (4.0±0.5±0.2±0.1)×10 −6 , compatible with the result for B(D 0 → K − π + µ + µ − ) [21] and with SM predictions [6].
We are grateful for the excellent luminosity and machine conditions provided by our PEP-II colleagues, and for the substantial dedicated effort from the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and kind hospitality. This work is supported by DOE and NSF (USA), NSERC (Canada), CEA and CNRS-IN2P3 I. Summary of signal candidate yields Nsig, average D 0 → K − π + e + e − reconstruction efficienciesˆ sig with statistical uncertainties, and B(D 0 → K − π + e + e − ) by m(e + e − ) mass region. The branching fraction uncertainties are statistical, systematic, and due to the known precision of B(D 0 → K − π + π + π − ), respectively. For m(e + e − ) ≤ 0.2 GeV/c 2 ,ˆ sig does not explicitly account for the photon in the decay π 0 → e + e − γ.