processes, such as angular distributions, charge-parity asymmetries or tests of lepton universality [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19].

Due to the universality of the electroweak interaction couplings to leptons in the SM, any differences in meson decay rates into final states with leptons of different generations arise purely from lepton mass effects. Consequently, comparing the branching fraction of a decay mode with muons in the final state to its counterpart with electrons provides a powerful probe of BSM effects, despite the dominance of long-distance contributions and large uncertainties in the SM prediction for the individual decay modes. Depending on the properties of possible BSM contributions in c → ul þ l -transitions, the validity of lepton universality could be broken at the percent level [6].

The

were first observed by the LHCb Collaboration [20], which subsequently studied their angular structure and charge-parity asymmetries [21,22]. Charge-conjugate decays are implied throughout the paper. The D 0 → h þ h -e þ e -decays offer the opportunity to probe the validity of lepton universality. However, electron final states are more challenging experimentally, in particular because the electrons lose a significant fraction of their energy as bremsstrahlung radiation while passing through the detector material. Until now, only upper limits on the branching fractions of D 0 → h þ h -e þ e -decays have been reported [23,24].

This paper reports the first search for the D 0 → π þ π -e þ e -and D 0 → K þ K -e þ e -decays at the LHCb experiment and opens the door for the first lepton universality tests with such decays. The analysis uses pp collision data collected between 2015 and 2018 at a center-of-mass energy ffiffi ffi s p ¼ 13 TeV, corresponding to an integrated luminosity of 6 fb -1 . The analysis method closely follows that of the LHCb analysis which observed the analogous decays with muons [20]; however the dielectron final state requires the development of a specific selection and reconstruction. The analysis uses D 0 mesons originating from D Ãþ → D 0 π þ decays with the D Ãþ mesons produced in the primary pp interaction.

The branching fractions of the decays of interest are measured in regions of the dielectron mass relative to the D 0 → K -π þ ½e þ e -ρ 0 =ω decay which has a branching fraction of ð4.0 AE 0.5Þ × 10 -6 in the dielectron-mass range 675-875 MeV=c 2 [25], where the notation indicates that the contribution from the ρ 0 =ω → e þ e -decay is dominant.

The LHCb detector is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5 and is described in detail in Refs. [26,27]. Events are selected by a trigger that consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage which applies a full event reconstruction [28]. The hardware trigger requires the presence of an energy deposit in the electromagnetic calorimeter that may or may not be compatible with originating from the signal candidate, or a muon signature with large transverse momentum which is compatible with originating from any particle in the event. A first stage of the software trigger selects events with either a charged particle that has both a significant transverse momentum and large impact parameter, defined as the minimum distance of the particle trajectory from any primary p p -collision vertex (PV), or with a two-track vertex satisfying a multivariate classifier based on geometric and kinematic criteria which identify the vertex as likely to originate from the decay of a long-lived heavy particle. In a second stage of the software trigger, candidate D 0 → h þ h ð0Þ-e þ e - decays are selected by combining four tracks that form a secondary vertex separated from any PV. All charged particles are required to have a significant impact parameter, as well as momentum p > 3 GeV=c and transverse momentum p T > 0.3 GeV=c. The D 0 candidate must have large transverse momentum, a reconstructed mass, mðD 0 Þ, in the range 1800-1950 MeV=c 2 (1700-2050 MeV=c 2 ) for data recorded in the years 2015-2016 (2017-2018) and its momentum vector must be aligned with the vector connecting the PVand the D 0 decay vertex. The larger mass range in the second period is due to a change in the software trigger that allows for improved background studies using information from the sidebands. When more than one PV is reconstructed, the one with respect to which the D 0 candidate has the lowest impact-parameter χ 2 is chosen, defined as the difference in the vertex-fit χ 2 reconstructed with and without the candidate. Finally, D Ãþ → D 0 π þ candidates are selected by combining the D 0 meson with a charged particle from the same PV that has p T > 120 MeV=c.

Subsequently, signal candidates are further selected offline by tightening the kinematic and geometric criteria applied in the trigger. A dedicated algorithm associates reconstructed bremsstrahlung photons to tracks identified as electrons; when a given photon is associated with both electron tracks, it is attached to one chosen randomly. Throughout the analysis, signal candidates are split into two periods with different data-taking conditions (2015-2016 and 2017-2018) and further divided into two categories: candidates in which neither electron has an associated bremsstrahlung cluster, and all other candidates, referred to as the no-brem and with-brem categories, respectively. A multivariate classifier is used to remove D 0 candidates containing one or more fake tracks. Stringent particleidentification criteria are then applied to all charged particles to suppress both combinatorial background, from unrelated charged particles, and cross feed backgrounds in which one type of D 0 → h þ h ð0Þ-e þ e -decay is misidentified as another. The vertex formed by the D 0 and π þ mesons is constrained to coincide with the PV and the momenta of the particles in the decay chain are updated accordingly. Only candidates with a difference between the reconstructed D Ãþ and D 0 masses, Δm, in the range 144-147 MeV=c 2 are considered. Further reduction of combinatorial background is achieved using a multivariate selection based on a boosted decision tree (BDT) [29,30]. The following features are used to discriminate signal from background: the momentum, transverse momentum, and impact parameter of the pion from the D Ãþ decay; the fit quality of the D 0 vertex and its separation from the PV; the angle between the D 0 momentum vector and the vector connecting the PV and the D 0 decay vertex; and the fit quality of the D Ãþ vertex. The BDT classifier is trained separately for

decays, using simulation samples [31,32] as a proxy for the signal and data candidates with mðD 0 Þ greater than 1900 MeV=c 2 as a proxy for the background.

Hadronic

where two pions are misidentified as electrons, constitute a major source of background, which is reduced by a multivariate electron-identification discriminant that combines information from the Cherenkov detectors, the calorimeters and the muon chambers. Finally, selection criteria based on the BDT response and on the electron identification discriminant are optimized by maximizing the figure-of-merit

is the signal efficiency, and N bkg is the total background yield in the mðD 0 Þ range 1700-1900 MeV=c 2 (1700-2050 MeV=c 2 ) for the no-brem (with-brem) category. The no-brem category has a significantly larger contamination from backgrounds in which hadrons are wrongly identified as electrons and therefore significantly worse sensitivity to the signals of interest. To facilitate an extrapolation of the yields of misidentified background from the no-brem to the with-brem category, the optimal selection thresholds obtained for the with-brem category are used for the no-brem category. Candidate D 0 → K -π þ ½e þ e -ρ 0 =ω decays are selected using the response of the BDT classifier trained on the

) signal when being used to normalize the branching fraction of

Only one candidate is kept at random if an event contains several signal or normalization candidate decays after the final selection, which happens in less than 0.5% of selected events. To avoid potential biases on the measured signal candidate yields, candidate decays in the mðD 0 Þ signal region were examined only after the analysis procedure was finalized.

Signal and normalization yields are measured with unbinned extended maximum-likelihood fits to the mðD 0 Þ distributions, in regions of the dielectron mass, mðe þ e -Þ, which is computed constraining the mass of the D 0 meson to its known value [34]. The fit results are shown for

candidates in Figs. 1 and 2, respectively. The fit results for D 0 → K -π þ ½e þ e -ρ 0 =ω can be found in the Supplemental Material of this paper [35]. For the signal modes, the mðe þ e -Þ regions are defined in Table I according to the presence of known intermediate resonances and aligned with Ref. [21]. Thus, the lowest mðe þ e -Þ region starts at 2 times the muon mass, m μ . All fits include four main components: signal, combinatorial background, peaking background from misidentified hadronic decays and background from partially reconstructed D 0 -meson decays, where at least one charged or neutral particle has not been reconstructed. The fit also accounts for the small fraction of

decays in which the dielectron pair is produced in one mðe þ e -Þ region and reconstructed in another. All fits are validated to return unbiased results using large numbers of pseudoexperiments. The signal is described with a Crystal Ball [36] distribution for the nobrem category and with a double-sided Crystal Ball distribution for the with-brem category. In both cases, the distribution parameters are determined from simulation. The mass shape of the peaking background is determined using a Bukin distribution [37] fitted to simulated samples of

, where the D 0 mass is calculated by assigning both pion candidates the electron mass. The peaking background in the with-brem category has a small, but not negligible, yield and a shape which is very similar to the signal. Its yield in the with-brem category is extrapolated from the no-brem category using

, η, ρ 0 =ω, ϕ and high-mðe þ e -Þ regions in the (left, 0γ) no-brem and (right, ≥ 1γ) with-brem categories. Fit projections are also shown. the relevant efficiencies from simulation. The combinatorial background is described by a first-order polynomial function with its slope determined from data candidates with Δm > 150 MeV=c 2 and mðD 0 Þ > 1900 MeV=c 2 . All shape parameters of the signal, peaking background and combinatorial background are fixed. Partially reconstructed backgrounds are modeled using a Bukin distribution fitted to simulated events in the low mðe þ e -Þ region (< 525 MeV=c 2 ), and an exponential function with its shape parameter determined in the fit to the selected candidates in other mðe þ e -Þ ranges. Alternative parametrizations are studied as a source of systematic uncertainty. The yields of each component are allowed to vary in the fits, which are performed simultaneously between the bremsstrahlung categories, periods of data-taking, and dielectron-mass regions. The yield of the combinatorial background is constrained to the yield determined in the Δm sideband 155-165 MeV=c 2 extrapolated to the signal region.

The resulting signal yields and significances with respect to zero, including statistical and systematic uncertainties, are reported in Table I. Significances exceeding 5 standard deviations are reported for D 0 → π þ π -e þ e -decays, where the two electrons are consistent with coming from an intermediate ϕ or ρ 0 =ω meson.

The signal yields,

region i are used to compute the branching fractions as

where

is the yield of the normalization mode, which is determined to be 820 AE 39 (875 AE 40) after applying the selection optimized for

corresponds to the ratio of geometrical acceptances, reconstruction and selection efficiencies of the signal relative to the normalization decays.

The efficiencies are determined using simulated events that are corrected to account for known differences between data and simulation. A particular challenge is the unknown amplitude composition of the decays under study. The signal decays are simulated with an incoherent sum of resonant and nonresonant dimuon and dihadron components. Samples of background-subtracted and efficiency-corrected D 0 → h þ h ð0Þ-μ þ μ -decays are used to correct the five-dimensional decay model of the corresponding dielectron mode. In addition, particleidentification, hardware-trigger, and tracking efficiencies as well as the dielectron-mass resolution are corrected using dedicated control channels in data.

Systematic uncertainties related to the determination of the yields arise due to limited knowledge of the various fit components and are evaluated using pseudoexperiments, where alternative fit models are tested. These include variations of the signal and background shape parameters within the uncertainties determined from the fits to simulation, and a signal shape obtained with a modified electron momentum resolution. In addition, a Bukin distribution is tested as an alternative distribution to represent the signal. Since not all selection criteria can be applied to background simulation because of limited simulated sample sizes, the misidentified background shapes are recomputed using an alternative set of selection criteria. The combinatorial background shape is determined in alternative Δm ranges. The fraction of signal decays migrating into other regions of dielectron mass is varied within its uncertainty. Furthermore, fit components are reevaluated using simulated samples corresponding to different datataking years. The dominant systematic uncertainty arising from assumptions in the fit is determined by fully neglecting all partially reconstructed backgrounds. The impact of neglecting further misidentified backgrounds is found to be negligible everywhere except in the signal ρ 0 =ω dielectronmass region, where a contribution from D 0 → K -π þ e þ e - decays misreconstructed as the signal are found, and the appropriate systematic uncertainty is computed. The statistical uncertainty on the normalization yield leads to a relative systematic uncertainty of 4.8% (4.5%) for

Systematic uncertainties affecting the efficiency ratio include residual data-simulation differences and limitations in the data-driven methods used to determine the particle-identification, tracking and trigger efficiencies. Uncertainties are evaluated directly on the efficiency ratio to take into account cancellations caused by similarities between the signal and normalization modes. The uncertainty due to finite simulated sample sizes is evaluated using a bootstrapping technique [38]. The impact of different detector occupancies in data and simulation is evaluated by recomputing the efficiencies using an additional correction accounting for the deviations. The efficiency ratios are recomputed varying the parameters used to

TABLE

signal decays and their significance, S, in units of Gaussian standard deviations, with respect to zero.

- Low mass 2m μ -525 37 AE 13 2.8σ η 525-565 10 AE 7 1 .6σ ρ 0 =ω 565-950 97 AE 21 5.5σ ϕ 950-1100 100 AE 18 8.1σ High mass > 1100 30 AE 11 2.9σ

- Low mass 2m μ -525 4 AE 8 1 .2σ η 525-565 1 AE 2 1 .1σ ρ 0 =ω > 565 12 AE 7 2 .2σ

smear the dielectron mass within their statistical uncertainties, and varying the resonant models of the signal decays within their uncertainties. In addition, the efficiency ratios are recomputed after artificially injecting an additional nonresonant component, representing unknown shortdistance physics, into the simulated signal distributions at 10% of the total number of events. Relative variations of the efficiency ratio between approximately 7% and 14% are found in the dielectron-mass regions and taken as systematic uncertainties. The impact of charm hadrons produced in beauty hadron decays is evaluated and found to be negligible.

To summarize, the biggest systematic uncertainties arise from finite simulated sample sizes, the limited knowledge of the signal-decay resonant structure, and the models used in the fit. The total systematic uncertainty depends on the decay mode and dielectron-mass region and ranges between approximately 10% and 70% of the corresponding statistical precision [35].

Measured branching fractions for II if the statistical significance exceeds 3 standard deviations, where the first uncertainty accounts for the statistical component, the second for the systematic, and the third corresponds to a 13.7% relative uncertainty on [25]. The Feldman-Cousins approach [39] is used to report the statistical uncertainties on the branching fractions. Upper limits are derived in all other mðe þ e -Þ regions using the CL s method including the effects due to the systematic uncertainties [40,41]. For comparison, the known branching fractions of [20], are also shown. Integrating over the dielectron-mass regions defined in Table II, and accounting for correlations [35], the total branching fraction for D 0 → π þ π -e þ e -decays with mðe þ e -Þ greater than 2 times the muon mass is measured to be

where the uncertainties are statistical, systematic and due to the limited knowledge of the normalization-mode branching fraction, respectively. The result is consistent with the SM expectations [5,6] and with the branching fraction of D 0 → π þ π -μ þ μ -of ð9.6 AE 1.2Þ × 10 -7 [20]. No total branching fraction is quoted for D 0 → K þ K -e þ e - decays as no significant signal is observed in any of the dielectron-mass regions.

In summary, a study of the

decays is performed in regions of the dielectron mass using p p collisions collected by the LHCb experiment at ffiffi ffi s p ¼ 13 TeV. The decay

is observed for the first time when requiring that the two electrons are consistent with coming from the decay of a ϕ or ρ 0 =ω mesons. No evidence is found for the

decay and world-best limits are set on its branching fraction which improves previous best limits by 2 orders of magnitude. The results are compared to, and found to be consistent with, the branching fractions of the

decays recently measured by LHCb, and confirm lepton universality at the current precision.

We express our gratitude to our colleagues in the CERN accelerator departments for the excellent performance of the LHC. We thank the technical and administrative staff at the LHCb institutes. We acknowledge support from CERN and from the national agencies: ARC (Australia); CAPES, CNPq, FAPERJ and FINEP (Brazil); MOST and NSFC (China); CNRS/IN2P3 (France); BMBF, DFG and MPG TABLE II. Branching fractions of (top) D 0 → π þ π -e þ e -and (bottom) D 0 → K þ K -e þ e -decays in different ranges of dielectron mass. Shown in comparison on the right are the branching fractions of

decays in the corresponding dimuon-mass ranges from Ref. [20]. The uncertainties are statistical, systematic and due to the limited knowledge of the normalization branching fraction. The reported upper limits correspond to 90% (95%) confidence level. The correlations between the various dielectron-mass ranges are reported in the Supplemental Material [35].

20] Low mass 2m μ -525 <4.8ð5.4Þ 0.78 AE 0.19 AE 0.05 AE 0.08 η 525-565 <2.3ð2.7Þ <0.24ð0.28Þ ρ 0 =ω 565-950 4.5 AE 1.0 AE 0.7 AE 0.6 4 .06 AE 0.33 AE 0.21 AE 0.41 ϕ 950-1100 3.8 AE 0.7 AE 0.4 AE 0.5 4 .54 AE 0.29 AE 0.25 AE 0.45 High mass >1100 <2.0ð2.2Þ <0.28ð0.33Þ D

20] Low mass 2m μ -525 <1.0ð1.1Þ 0.26 AE 0.12 AE 0.02 AE 0.03 η 525-565 <0.4ð0.5Þ <0.07ð0.08Þ ρ 0 =ω >565 <2.2ð2.5Þ 1.20 AE 0.23 AE 0.07 AE 0.12