Search for the rare decays $D\to h(h')e^+e^-$

We search for rare decays of $D$ mesons to hadrons accompany with an electron-positron pair (h(h')$e^+e^-$), using an $e^+e^-$ collision sample corresponding to an integrated luminosity of 2.93 fb$^{-1}$ collected with the BESIII detector at $\sqrt{s}$ = 3.773 GeV. No significant signals are observed, and the corresponding upper limits on the branching fractions at the $90\%$ confidence level are determined. The sensitivities of the results are at the level of $10^{-5} \sim 10^{-6}$, providing a large improvement over previous searches.


I. INTRODUCTION
In the standard model (SM), the decay of a D meson into hadrons accompanied by a lepton pair proceeds via the quark process c → ul þ l − (l ¼ e or μ). This is known as a flavor changing neutral current (FCNC) process, which is forbidden at tree level in the SM. It can happen only through a loop diagram because of the suppression of the Glashow-Iliopoulos-Maiani (GIM) mechanism [1], leading to a very small branching fraction (BF) theoretically, which would not exceed the level of 10 −9 [2][3][4]. Compared to similar FCNC processes in Band K-meson decays, the GIM suppression in FCNC decays of the D meson is much stronger, as better diagram cancellation occurs due to the down-type quarks involved. However, possible new physics (NP) beyond the SM can significantly increase the decay rates of these short distance (SD) processes. Hence, they can serve as clean channels in experiments to search for NP [2,3].
However, these D-meson-decay rates are also contributed by long distance (LD) effects through (virtual) vector meson (V ð⋆Þ ) decays, like D → hV ð⋆Þ ; V ð⋆Þ → l þ l − , even above the level of 10 −6 [3,4]. Therefore, FCNC processes are potentially overshadowed by LD effects. In such case, a measurement of the angular dependence or CP asymmetry is required to figure out the SD effects and to test the SM prediction.
In recent years, the four-body decays of D 0 mesons with a μ þ μ − pair in final state, i.e., D 0 → K − π þ μ þ μ − , K − K þ μ þ μ − and π − π þ μ þ μ − , have been observed at LHCb [5,6], with the decay rates at the level of 10 −7 , indicating significant LD contributions. However, no evidence for the e þ e − modes has yet been reported. The three-body and four-body decays of D 0 mesons involving e þ e − pairs were searched for by the CLEO and E791 Collaborations [7,8]. The current upper limits (ULs) on their branching fractions at the 90% confidence level (CL) are at the level of 10 −4 -10 −5 . The analogous D þ decays are less well studied, and only three-body decays have been searched for by the BABAR and LHCb Collaborations [9].
In this paper, using an e þ e − collision sample corresponding to an integrated luminosity of 2.93 fb −1 [10] collected with the BESIII detector at ffiffi ffi s p ¼ 3.773 GeV, we perform a search for the rare decays of D → hðh ð0Þ Þe þ e − , where h ð0Þ are hadrons. To avoid possible bias, a blind analysis is carried out based on Monte Carlo (MC) simulations to validate the analysis strategy, the results are opened only after the analysis strategy is fixed.

II. THE BESIII DETECTOR AND MC SIMULATION
The Beijing Spectrometer (BESIII) detects e þ e − collisions in the double-ring collider BEPCII. BESIII is a general-purpose detector [11] with 93% coverage of the full solid angle. From the interaction point (IP) to the outside, BESIII is equipped with a main drift chamber (MDC) consisting of 43 layers of drift cells, a time-of-flight (TOF) counter with double-layer scintillator in the barrel part and single-layer scintillator in the end-cap part, an electromagnetic calorimeter (EMC) composed of 6240 CsI (Tl) crystals, a superconducting solenoid magnet providing a magnetic field of 1.0 T along the beam direction, and a muon counter containing multilayer resistive plate chambers installed in the steel flux-return yoke of the magnet. The MDC spatial resolution is 135 μm and the momentum resolution is 0.5% for a charged track with transverse momentum of 1 GeV=c. The energy resolution for a photon at 1 GeV in the EMC is 2.5% in the barrel region and 5.0% in the endcap region. More details of the spectrometer can be found in Ref. [11].
Monte Carlo simulation serves to estimate the detection efficiencies and to understand background contamination. High statistics MC samples are generated with a GEANT4based [12] software package, which includes the descriptions of the geometry of the spectrometer and interactions of particles with the detector materials. KKMC [13] is used to model the beam energy spread and the initial-state radiation (ISR) in the e þ e − annihilations. The "inclusive" MC samples consist of the production of DD pairs with quantum coherence for all neutral D modes, the non-DD decays of ψð3770Þ, the ISR production of low mass ψ states, and continuum processes. Known decays recorded by the Particle Data Group (PDG) [9] are simulated with EVTGEN [14] and the unknown decays with LUNDCHARM model [15]. The final-state radiation (FSR) of charged tracks is taken into account with the PHOTOS package [16]. The equivalent luminosity of the inclusive MC samples is about 10 times that of the data. The signal processes are generated using the phase space model (PHSP) of EVTGEN. For each signal channel, 200000 events are simulated.

III. DATA ANALYSIS
Since the center-of-mass energy of 3.773 GeV is close to the DD mass threshold, the pair of D þ D − or D 0D0 mesons is produced nearly at rest without any additional hadrons. Hence, it is straightforward to use a double tagging approach [17] to measure absolute BFs, based on the following equation Here, i denotes the different single-tag (ST) modes of hadronic decays, and n i tag is the yield of theD meson of ST tag mode i. n sig;tag is the number of D rare decay candidate events in which a STD meson is detected, so called doubletag (DT) events. Finally, ε i tag and ε i sig;tag are the corresponding ST and DT detection efficiencies. The average signal efficiency over different ST modes can be calculated to be where n tag is the total number of ST events n tag ¼ P i n i tag . Note that in this paper, charge conjugated modes are always implied.
Momenta and impact parameters of charged tracks are measured by the MDC. Charged tracks (except for those of K 0 S decays) are required to satisfy jcos θj < 0.93, where θ is the polar angle with respect to the beam axis, and have a distance of closest approach to the interaction point (IP) within AE10 cm along the beam direction and within AE1 cm in the plane perpendicular to the beam axis. Particle identification (PID) is implemented by combining the specific energy loss (dE=dx) in the MDC and the time of flight measured from the TOF to form PID confidence levels (CL) for each particle hypothesis. For a charged πðKÞ candidate, the CL of the πðKÞ hypothesis is required to be larger than that of the KðπÞ hypothesis.
Photon candidates are reconstructed from clusters of deposited energy in the EMC. The energies of photon candidates must be larger than 25 MeV for jcos θj < 0.8 (barrel) or 50 MeV for 0.86 <jcos θj < 0.92 (end cap). To suppress fake photons due to electronic noise or beam background, the shower time must be less than 700 ns from the event start time [18]. The photon candidates are required to be at least 20°away from any charged track.
The π 0 candidates are reconstructed from pairs of photons of which at least one is reconstructed in the barrel. The invariant mass of the photon pair, M γγ , is required to lie in the range ð0.115; 0.150Þ GeV=c 2 . To improve the resolution, we further constrain the invariant mass of each photon pair to the nominal π 0 mass by a kinematic fit, and the updated four-momentum of the π 0 will be used in the further analysis.
The K 0 S candidates are reconstructed via K 0 S → π þ π − using a vertex-constrained fit to all pairs of oppositely charged tracks, without PID requirements. The distance of closest approach of a charged track to the IP is required to be less than AE20 cm along the beam direction, without any requirement in the transverse plane. The χ 2 of the vertex fit is required to be less than 100. The invariant mass of the π þ π − pair, M π þ π − , is required to be within ð0.487; 0.511Þ GeV=c 2 , corresponding to three times the experimental mass resolution.
Two variables, the beam-constrained mass, M tag BC , and the energy difference, ΔE tag , which are defined as are used to identify the tag candidates. Here, ⃗pD and ED are the momentum and energy of the STD candidate in the rest frame of the initial e þ e − system, and E beam is the beam energy. Signal events peak around the nominalD mass in the M tag BC distribution and around zero in the ΔE tag distribution. The boundaries of the ΔE tag requirements are determined from MC simulation, and set at approximately (−3σ, 3σ) for the modes with only charged tracks in final state, and (−4σ; þ3.5σ) for those including a π 0 in the final state, due to the asymmetric ΔE distribution. Here, σ is the standard deviation of ΔE tag . In each event, only the combination with the smallest jΔE tag j is kept for each ST mode.
After applying the ΔE tag requirements as listed in Table I for the different ST modes, the M tag BC distributions are shown in Fig. 1. The corresponding ST yields are extracted by performing maximum likelihood fits to the M tag BC distribution, where in each mode the signal is modeled with a MCderived signal shape convolved with a smearing Gaussian function representing the resolution difference between data and MC simulation, and the backgrounds are modeled with an ARGUS function [19]. Based on the fit results, the total ST yields found in data are summarized in Table I together with the MC-determined detection efficiencies.
Events with a ST candidate fulfilling the additional requirement of M tag BC to be within ð1.863; 1.879Þ GeV=c 2 for charged D modes and ð1.858; 1.874Þ GeV=c 2 for neutral D modes are used to search for signal candidates as described in the following.

B. Signal event selection and yields
Signal candidates of D þ decaying to π þ π 0 e þ e − , K þ π 0 e þ e − , K 0 S π þ e þ e − , and K 0 S K þ e þ e − , and D 0 decaying to K − K þ e þ e − , π þ π − e þ e − , K − π þ e þ e − , π 0 e þ e − , ηe þ e − , ωe þ e − , and K 0 S e þ e − are searched for in the remaining charged tracks and showers recoiling against the STD mesons. The selection criteria for the charged tracks and neutral showers are the same as those used in the ST event selection. Positrons and electrons are distinguished from other charged particles by combining the dE=dx, TOF and EMC information. The determined particle identification CL, L, is required to satisfy LðeÞ > 0 and LðeÞ=ðLðeÞ þ LðπÞ þ LðKÞÞ > 0.8. Furthermore, the energy deposited in the EMC divided by the momentum measured in the MDC, E=pc, is required to be larger than 0.8 for either the positron or electron. By studying the inclusive MC samples, we find that the selected e þ e − pairs dominantly originate from γ-conversion events, where the photons are from the decay of intermediate states. To suppress these backgrounds, the vertex of the e þ e − pair is reconstructed [20,21], and the distance from the IP to the reconstructed vertex in the x − y plane R xy is required to be out of range (2.0, 8.0) cm, where the γ-conversion occurs.
To veto the contribution from D → hðh ð0Þ Þϕ; ϕ → e þ e − , the e þ e − invariant mass M e þ e − is required to be outside of the ϕ mass region, defined as ð0.935; 1.053Þ GeV=c 2 .
An η candidate is reconstructed via its γγ decay mode by requiring M γγ within ð0.505; 0.570Þ GeV=c 2 . A kinematic fit constraining M γγ to the nominal η mass is applied, and the candidate with the smallest χ 2 is kept under the requirement χ 2 < 20. Similarly, candidate π 0 decaying into γγ are selected by requiring M γγ within ð0.110; 0.155Þ GeV=c 2 . A kinematic fit constraining M γγ to the nominal π 0 mass is performed. The candidate with the smallest χ 2 is kept and is required to satisfy χ 2 < 20. An ω candidate is reconstructed with its π þ π − π 0 decay mode, by requiring the three-pion invariant mass M π þ π − π 0 to be within ð0.720; 0.840Þ GeV=c 2 . For the K 0 S candidates, in addition to the same criteria as used in ST event selection, we further require L=σ L > 2, where L is the measured K 0 S flight distance and σ L is the corresponding uncertainty.
Similar to the ST selection, ΔE and M BC for the signal candidates of the rare D decays in DT events, denoted as ΔE sig and M sig BC , are calculated. For each signal mode, ΔE sig   is required to be within 3σ of the nominal value, as listed in Table II, and only the combination with the smallest jΔE sig j is kept. The M sig BC distributions of the surviving events are shown in Fig. 2, where no significant excess over the expected backgrounds is observed. The number of remaining signal candidates, n obs , is counted in the M sig BC signal regions and listed in Table II. The corresponding DT detection efficiencies and the average signal efficiencies ε sig over different ST modes are given in Table III. The BFs of the rare decays will be determined by subtracting the background contributions.
The backgrounds are separated into two categories: events with a wrong ST candidate, and events with a correct ST but wrong signal candidate, which dominantly originate from the γ-conversion process. The former background can be estimated with the surviving events in the ST sideband (SB) region of M tag BC distribution, which is defined as ð1.830; 1.855Þ GeV=c 2 forD 0 decays and ð1.830; 1.860Þ GeV=c 2 for D − decays. The corresponding number of wrong-ST background events, n bkg1 , is estimated with the number of events in the SB region (n SB bkg1 ) normalized by a scale factor f, which is the ratio of the integrated numbers of background events in the signal and SB regions. The scale factor f is found to be 0.466 AE 0.001 for the D þ decays and 0.611 AE 0.001 for the D 0 decays, respectively, where the uncertainty is statistical only. The wrong-ST background is expected to follow a Poisson (P) distribution with central value of n bkg1 · f. The background from misreconstructed signal is estimated with the D þ D − and D 0D0 events in the inclusive MC samples by subtracting the wrong STevents, and the corresponding number of events is expected to follow a Gaussian distribution (G), with central value n MC bkg2 and standard deviation σ MC bkg2 . The relevant numbers are summarized in Table II.

IV. SYSTEMATIC UNCERTAINTIES
With the DT technique, the systematic uncertainties in the BF measurements due to the detection and reconstruction of the STD mesons mostly cancel, as shown in Eq. (1). For the signal side, the following sources of systematic uncertainties, as summarized in Table IV, are considered. All of these contributions are added in quadrature to obtain the total systematic uncertainties.
The uncertainties of tracking and PID efficiencies for K AE and π AE are studied with control samples of DD favored hadronic modes [22]. We assign an uncertainty of 1.0% per track for the tracking and 0.5% for the PID uncertainties. The tracking and PID efficiency for e AE detection is studied using radiative Bhabha events, and the corresponding systematic uncertainty is evaluated by weighting according to the cos θ and transverse momentum distributions of the e AE tracks. The uncertainties for π 0 , η and K reconstructions are studied with control samples of DD events. An uncertainty of 2.0% is assigned for each π 0 , 1.5% for K 0 s , and 1.2% for η. The γ-conversion background is suppressed by a requirement on the distance from the reconstructed vertex of the e þ e − pair to the IP. The uncertainty due to this requirement is studied using a sample of J=ψ → π þ π − π 0 with π 0 → γe þ e − [21]. The relative difference of the efficiency between data and MC simulation is 1.8%, and is assigned as the uncertainty.
The estimated signal detection efficiencies depend on the MC simulations, which is assumed to be distributed uniformly in momentum phase space. However, theoretically there are nontrivial contributions from LD effects [4]. Alternative MC samples with LD models, in which the e þ e − pairs originate from vector mesons are also generated to estimate the signal detection efficiencies. The resultant changes on the detection efficiencies are assigned as the systematic uncertainty. The BF uncertainty for the intermediate states decays of the neutral mesons, B inter , are assigned according to the world average values [9].

V. THE UPPER LIMITS ON BRANCHING FRACTIONS
To calculate the ULs on the BFs for the signal decays, we use a maximum likelihood estimator, extended from the profile likelihood method [23]. For the detection efficiency, we assume it follows a Gaussian distribution, whose mean and width are MC-determined efficiency ε MC sig and its absolute systematic uncertainty ε MC sig · σ MC ε , respectively, and σ MC ε includes the relative statistical and systematic uncertainties as given in Table III and Table IV. So the joint likelihood is     L ¼ Pðn obs ; n tag · B · ε sig þ n bkg1 þ n bkg2 Þ · Gðε sig ; ε MC sig ; ε MC sig · σ MC ε Þ · Pðn SB bkg1 ; n bkg1 · fÞ · Gðn bkg2 ; n MC bkg2 ; σ MC bkg2 Þ: ð2Þ Based on the Bayesian method, we use the likelihood distribution as a function of the signal BF B, with variations of the other parameters n bkg1 , n bkg2 , and ε sig , as the probability function. Note that the ST yields, n tag , are taken as the truth ones, as their uncertainties are negligible. The resultant likelihood distributions for all the signal modes are shown in Fig. 3, and the ULs on the signal BFs at the 90% CL are estimated by integrating the likelihood curves in the physical region of B ≥ 0. For D 0 → K − π þ e þ e − , the BF is determined to be ð2.5 AE 1.1Þ × 10 −5 with a significance of 2.6σ, where the uncertainty includes the statistical and systematic ones. Reference [4] predicts the BF of D 0 → K − π þ e þ e − , which is dominated by the LD bremsstrahlung and (virtual) resonance-decay contributions in the lower and upper regions, respectively, to exceed 0.99 × 10 −5 in the lower M e þ e − region, adding up to 1.6 × 10 −5 in the whole region. Therefore, we divide the M e þ e − distribution into three regions and determine the BFs in the individual regions. All these results are listed in Table V, and are all within the SM predictions.  signals are observed, and the corresponding ULs on the decay rates are determined at the 90% CL, as shown in Table V. For the four-body D þ decays, the searches are performed for the first time. The reported ULs of the D 0 decays are improved in general by a factor of 10, compared to previous measurements [9]. All the measured ULs on the BFs are above the SM predictions [3,4], which include both LD and SD contributions.