Observation of $D^0$ meson decays to $\pi^+\pi^-\mu^+\mu^-$ and $K^+K^-\mu^+\mu^-$ final states

The first observation of the $D^0 \to \pi^+\pi^-\mu^+\mu^-$ and $D^0 \to K^+K^-\mu^+\mu^-$ decays is reported using a sample of proton-proton collisions collected by LHCb at a center-of-mass energy of 8$\,$TeV, and corresponding to 2$\,$fb$^{-1}$ of integrated luminosity. The corresponding branching fractions are measured using as normalization the decay $D^0 \to K^- \pi^+[\mu^+\mu^-]_{\rho^0/\omega}$, where the two muons are consistent with coming from the decay of a $\rho^0$ or $\omega$ meson. The results are $\mathcal{B}(D^0 \to \pi^+\pi^-\mu^+\mu^-)=(9.64\pm0.48\pm0.51\pm0.97)\times10^{-7}$ and $\mathcal{B}(D^0 \to K^+K^-\mu^+\mu^-)=( 1.54\pm0.27\pm0.09\pm0.16)\times10^{-7}$, where the uncertainties are statistical, systematic, and due to the limited knowledge of the normalization branching fraction. The dependence of the branching fraction on the dimuon mass is also investigated.

Decays of charm hadrons into final states containing dimuon pairs may proceed via the short-distance c → uμ þ μ − flavor-changing neutral-current process, which in the standard model can only occur through electroweakloop amplitudes that are highly suppressed by the Glashow-Iliopoulos-Maiani mechanism [1]. If dominated by these short-distance contributions, the inclusive D → Xμ þ μ − branching fraction, where X represents one or more hadrons, is predicted to be Oð10 −9 Þ [2] and can be greatly enhanced by the presence of new particles, making these decays interesting for searches for physics beyond the standard model. However, long-distance contributions occur through tree-level amplitudes involving intermediate resonances, such as D → XVð→ μ þ μ − Þ, where V represents a ρ 0 , ω or ϕ vector meson, and can increase the standard model branching fraction up to Oð10 −6 Þ [2-4]. The sensitivity to the short-distance amplitudes is greatest for dimuon masses away from resonances, though resonances populate the entire dimuon-mass spectrum due to their long tails. Additional discrimination between shortand long-distance contributions can be gained by studying angular distributions and charge-parity-conjugation asymmetries, which in scenarios beyond the standard model could be as large as Oð1%Þ [4][5][6][7][8][9]. Decays of D 0 mesons to four-body final states ( Fig. 1) are particularly interesting in this respect as they give access to a variety of angular distributions. These decays were searched for by the Fermilab E791 Collaboration and upper limits were set on the branching fractions in the range 10 −5 − 10 −4 at the 90% confidence level (C.L.) [10]. More recently, a search for nonresonant D 0 → π þ π − μ þ μ − decays (the inclusion of charge-conjugate decays is implied) was performed by the LHCb Collaboration using 7 TeV pp-collision data corresponding to 1 fb −1 of integrated luminosity [11]. An upper limit of 5.5 × 10 −7 at the 90% C.L. was set on the branching fraction due to short-distance contributions, assuming a phase-space decay.
This Letter reports the first observation of D 0 → π þ π − μ þ μ − and D 0 → K þ K − μ þ μ − decays using data collected by the LHCb experiment in 2012 at a center-of-mass energy ffiffi ffi s p ¼ 8 TeV and corresponding to an integrated luminosity of 2 fb −1 . The analysis is performed using D 0 mesons originating from D Ãþ → D 0 π þ decays, with the D Ãþ meson produced directly at the primary pp-collision vertex (PV). The small phase space available in this decay allows for a large background rejection, which compensates for the reduction in signal yield compared to inclusively produced D 0 mesons. The signal is studied in regions of dimuon mass, mðμ þ μ − Þ, defined according to the known resonances. For D 0 → π þ π − μ þ μ − decays these regions are (low-mass) < 525 MeV=c 2 , (η) 525-565 MeV=c 2 , (ρ 0 =ω) 565-950 MeV=c 2 , (ϕ) 950-1100 MeV=c 2 , and (highmass) > 1100 MeV=c 2 . The same regions are considered for D 0 → K þ K − μ þ μ − decays, with the exception of the ϕ and high-mass regions, which are not present because of the reduced phase space, and the ρ 0 =ω region, which extends from 565 MeV=c 2 up to the kinematic limit. In the regions where a signal is observed a measurement of the branching fraction is provided, otherwise 90% and 95% C.L. upper limits are set; no attempt is made to distinguish between the short-and long-distance contributions in each dimuon-mass region. The branching fraction is measured using as a normalization the D 0 → K − π þ ½μ þ μ − ρ 0 =ω decay in the dimuon-mass range 675-875 MeV=c 2 , where the contribution from the ρ 0 =ω → μ þ μ − decay is dominant. The D 0 → K − π þ ½μ þ μ − ρ 0 =ω branching fraction was recently measured to be ð4.17 AE 0.42Þ × 10 −6 [12] and provides a more precise normalization than that used in the previous LHCb search [11]. The LHCb detector is a single-arm forward spectrometer [13,14]. It includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp-interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. Particle identification is provided by two ring-imaging Cherenkov detectors, an electromagnetic and a hadronic calorimeter, and a muon system composed of alternating layers of iron and multiwire proportional chambers.
Events are selected online by a trigger that consists of a hardware stage, which is based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction [15]. The hardware trigger requires the presence in the event of a muon with transverse momentum, p T , exceeding 1.76 GeV=c. A first stage of the software trigger selects events with a charged particle of p T > 1.6 GeV=c and significant impact parameter, defined as the minimum distance of the particle trajectory from any PV, or alternatively with p T > 1 GeV=c if the particle has associated hits in the muon system. In a second stage of the software trigger, dedicated algo- where h is either a kaon or a pion, from combinations of four tracks, each having momentum p > 3 GeV=c and p T > 0.5 GeV=c, that form a secondary vertex separated from any PV. Two oppositely charged particles are required to leave hits in the muon system and the scalar sum of their p T is required to exceed 3 GeV=c. The mass of the D 0 candidate, mðD 0 Þ, has to be in the range 1800-1940 MeV=c 2 and its momentum must be aligned with the vector connecting the primary and secondary vertices.
In the offline analysis, D 0 candidates satisfying the trigger requirements are further selected through particleidentification criteria placed on their decay products. They are then combined with a charged particle originating from the same PV and having p T > 120 MeV=c, to form a D Ãþ → D 0 π þ candidate. When more than one PV is reconstructed, the one with respect to which the D 0 candidate has the lowest impact-parameter significance is chosen. The vertex formed by the D 0 and π þ mesons is constrained to coincide with the PV and the difference between the D Ãþ and D 0 masses, Δm, is required to be in the range 144.5-146.5 MeV=c 2 . A multivariate selection based on a boosted decision tree (BDT) [16,17] with gradient boosting [18] is then used to suppress background from combinations of unrelated charged particles. The features used by the BDT to discriminate signal from this combinatorial background are as follows: the momentum and transverse momentum of the pion from the D Ãþ decay, the smallest impact parameter of the D 0 decay products with respect to the PV, the angle between the D 0 momentum and the vector connecting the primary and secondary vertices, the quality of the secondary vertex, its separation from the PV, and its separation from any other track not forming the D Ãþ candidate. The BDT is trained separately for D 0 → π þ π − μ þ μ − and D 0 → K þ K − μ þ μ − decays, due to their different kinematic properties, using simulated [19,20] decays as signal and data candidates with mðD 0 Þ between 1890 and 1940 MeV=c 2 as background. To minimize biases on the background classification, the training samples are further randomly split into two disjoint subsamples. The classifier trained on one sample is applied to the other, and vice versa. Another source of background is due to the hadronic four-body decays D 0 → π þ π − π þ π − and D 0 → K þ K − π þ π − , where two pions are misidentified as muons. The misidentification occurs mainly when the pions decay in flight into a muon and an undetected neutrino. Although this process is relatively rare, the large branching fractions of the hadronic modes produce a peaking background which is partially suppressed by a multivariate muon-identification discriminant that combines the information from the Cherenkov detectors, the calorimeters and the muon chambers. Thresholds on the BDT response and on the muon-identification discriminant are optimized simultaneously by maximizing week ending 3 NOVEMBER 2017 181805-2 signal efficiency and N bkg is the sum of the expected combinatorial and peaking background yields in the mðD 0 Þ range 1830-1900 MeV=c 2 (signal region). Candidate D 0 → K − π þ ½μ þ μ − ρ 0 =ω decays are selected using the response of the BDT trained on the D 0 → π þ π − μ þ μ − signal, when they are used as normalization for the measurement of BðD 0 → π þ π − μ þ μ − Þ, and that of the BDT trained on the D 0 → K þ K − μ þ μ − signal, when used as normalization for BðD 0 → K þ K − μ þ μ − Þ. After selection, a few percent of the events contain multiple candidates, of which only one is randomly selected if they share at least one final-state particle. To avoid potential biases on the measured quantities, candidate decays in the mðD 0 Þ signal region are examined only after the analysis procedure has been finalized, with the exception of those populating the ρ 0 =ω and ϕ dimuon-mass regions of the are measured with unbinned extended maximum likelihood fits to the mðD 0 Þ distributions (Figs. 2 and 3, respectively). The fits include three components: signal, peaking background from misidentified hadronic decays, and combinatorial background. The signal is described with a Johnson's S U distribution [22] with parameters determined from simulation. To account for known differences between data and simulation, the means and widths of the signal distributions are corrected using scaling factors adjusted on the normalization channel. The mass shape of the peaking background is determined using separate data samples of D 0 → h þ h ð0Þ− π þ π − decays where the D 0 mass is calculated assigning the muon-mass hypothesis to two oppositely charged pions. The combinatorial background is described by an exponential function, which is determined from data candidates with Δm between 150 and 160 MeV=c 2 that fail the BDT selection. All shape parameters are fixed and only the yields are allowed to vary in the fits, which are performed separately in each mðμ þ μ − Þ range.
The resulting signal yields are reported in Table I. No fit is performed in the η region of the D 0 → K þ K − μ þ μ − dimuon-mass spectrum, where only two candidates are observed. An excess of candidates with respect to the background-only hypothesis is seen with a significance above three standard deviations in all dimuon-mass ranges with the exception of the η region of both decays and the high-mðμ þ μ − Þ region of D 0 → π þ π − μ þ μ − . The significances are determined from the change in likelihood from fits with and without the signal component.
The signal yields, where N K − π þ μ þ μ − is the yield of the normalization mode, which is determined to be 1971 AE 51 (1806 AE 48) after the selection optimized for The ratios of geometrical acceptances, and reconstruction and selection efficiencies of the signal relative to the normalization decays, Table I. They are determined using simulated events and corrected to account for known differences between data and simulation. In particular, particle-identification and hardwaretrigger efficiencies are measured from control channels in data.
Systematic uncertainties affect the determination of the signal and normalization yields, and of the efficiency ratio. For the determination of the yields, effects due to uncertainties on the mðD 0 Þ shapes are investigated. A possible dependence on the decay mode or on the mðμ þ μ − Þ range of the scaling factors, used to account for data-simulation differences, is quantified using fits to the D 0 → π þ π − ½μ þ μ − ϕ and D 0 → π þ π − ½μ þ μ − ρ 0 =ω data and is found to be negligible. To assess the impact of π → μν decays in flight, alternative shapes are tested for the D 0 → h þ h ð0Þ− π þ π − background by changing the muonidentification and the p T requirements on the misidentified pions. The largest observed variation in the ratio of D 0 → π þ π − ½μ þ μ − ϕ to D 0 → K − π þ ½μ þ μ − ρ 0 =ω yields (1.4%) is assigned as a systematic uncertainty for both h þ h − μ þ μ − modes and all dimuon-mass ranges. Changes in the shape of the peaking background introduced by the different trigger requirements used to select the hadronic decays are negligible. The fit to the data is repeated using alternative descriptions of the combinatorial background, determined from data sidebands defined by different BDT and Δm requirements, and results in negligible variations of the signal and normalization yields.
Systematic uncertainties affecting the efficiency ratio include data-simulation differences that are not accounted for and limitations in the data-driven methods used to determine the particle-identification and trigger efficiencies. The signal decays are simulated with an incoherent sum of resonant and nonresonant dimuon and dihadron components, while the resonant structure in data is unknown. A systematic uncertainty of 3.4% on the signal efficiency is determined by varying the relative fractions of these components. A systematic uncertainty of 1.0% on the efficiency ratio is assigned due to the criteria used in simulation to match the reconstructed and generated particles. Muon-and hadron-identification efficiencies are determined from data by weighting the kinematic properties of the calibration samples to match those of the signal samples. Variations of the choice of the binning scheme used in the weighting procedure change the efficiency ratio by up to 0.8%, which is taken as systematic uncertainty. The datadriven method that evaluates the hardware-trigger efficiency ratio is validated in simulation to be unbiased within 1.3%, which is assigned as a systematic uncertainty. The efficiencies of the BDT requirement for the simulated normalization and D 0 → π þ π − ½μ þ μ − ϕ decays are compared to those obtained from background-subtracted data. A difference in the efficiency ratio of 1.3% is observed and assigned as systematic uncertainty.
Finally, the statistical uncertainty on the normalization yield introduces a relative uncertainty of 2.6% (2.7%), which is propagated to the systematic uncertainty on the Table II reports the measured values and upper limits on the D 0 → π þ π − μ þ μ − and D 0 → K þ K − μ þ μ − branching fractions in the various ranges of mðμ þ μ − Þ, where the first uncertainty accounts for the statistical component, the second for the systematic, and the third corresponds to the 10% relative uncertainty on BðD 0 → K − π þ ½μ þ μ − ρ 0 =ω Þ [12]. The upper limits are derived using a frequentist approach based on a likelihood-ratio ordering method that includes the effects due to the systematic uncertainties [23,24]. For the η region of D 0 → K þ K − μ þ μ − , where no fit is performed, the limit is calculated assuming two signal TABLE I. Yields of (top) D 0 → π þ π − μ þ μ − and (bottom) D 0 → K þ K − μ þ μ − signal decays, their significance with respect to the background-only hypothesis, and ratio of efficiencies between signal and normalization decays (R i ϵ ) for each dimuon-mass region. The yield and the significance (S) are not reported for the η region of D 0 → K þ K − μ þ μ − , where only two candidates are observed.
candidates and zero background. Integrating over dimuon mass, and accounting for correlations [25], the total branching fractions are measured to be The two results have a correlation of 0.497 and are consistent with the standard model expectations [4]. In summary, a study of the D 0 → π þ π − μ þ μ − and D 0 → K þ K − μ þ μ − decays is performed in ranges of the dimuon mass using pp collisions collected by the LHCb experiment at ffiffi ffi s p ¼ 8 TeV. Significant signal yields are observed for the first time in several dimuon-mass ranges for both decays; the corresponding branching fractions are measured and found to be consistent with the standard model expectations [4]. For the dimuon-mass regions where no significant signal is observed, upper limits at 90% and 95% C.L. are set on the branching fraction. The total branching fractions are measured to be BðD where the uncertainties are statistical, systematic, and due to the limited knowledge of the normalization branching fraction. These are the rarest charm-hadron decays ever observed and are expected to provide better sensitivity to short-distance flavor-changing neutral-current contributions to these decays.
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: CAPES, CNPq, FAPERJ and FINEP (Brazil)    Branching fractions of (top) D 0 → π þ π − μ þ μ − and (bottom) D 0 → K þ K − μ þ μ − decays in different ranges of dimuon mass, where the uncertainties are statistical, systematic, and due to the limited knowledge of the normalization branching fraction. The reported upper limits correspond to 90% (95%) C.L. The correlations between the various dimuon-mass ranges are reported in the Supplemental Material [25].