Angular analysis of $D^0 \to \pi^+\pi^-\mu^+\mu^-$ and $D^0 \to K^+K^-\mu^+\mu^-$ decays and search for $CP$ violation

The first full angular analysis and an updated measurement of the decay-rate $CP$ asymmetry of the $D^0 \to \pi^+\pi^-\mu^+\mu^-$ and $D^0 \to K^+K^-\mu^+\mu^-$ decays are reported. The analysis uses proton-proton collision data collected with the LHCb detector at centre-of-mass energies of 7, 8 and 13 TeV. The data set corresponds to an integrated luminosity of 9 fb$^{-1}$. The full set of $CP$-averaged angular observables and their $CP$ asymmetries are measured as a function of the dimuon invariant mass. The results are consistent with expectations from the standard model and with $CP$ symmetry.

Rare charm decays with two oppositely charged leptons ( + − ) in the final state may proceed via the quark flavor-changing neutral-current (FCNC) process c → u + − and, as such, be sensitive to contributions from physics beyond the standard model (SM). They represent a unique probe to beyond-SM couplings to up-type quarks, which is complementary to recent studies of beauty quark b → s + − transitions, where a coherent pattern of deviations from the SM is emerging (e.g., see Refs. [1][2][3]). The loop-induced SM processes are more suppressed in charm than in the b-quark system due to the Glashow-Iliopoulos-Maiani mechanism [4]. The short-distance contributions to the inclusive D → Xµ + µ − branching fraction, where D denotes a neutral or charged D meson and X represents one or more hadrons, are predicted to be of O(10 −9 ) [5]. Sensitivity to FCNC processes via the measurement of branching fractions is limited due to the dominance of tree-level amplitudes involving intermediate resonances that subsequently decay into + − . These so-called long-distance contributions increase the SM branching fractions up to O(10 −6 ) [5][6][7][8]. Studies of angular distributions and charge-parity (CP ) asymmetries in the vicinity of intermediate resonances offer a access to observables with negligible theoretical uncertainties. These observables are sensitive to beyond-SM physics through the interference between long-and short-distance amplitudes. The values of these observables are negligibly small in the SM, but can reach the percent level in scenarios beyond the SM [7][8][9][10][11][12][13][14][15][16][17][18].
The LHCb Collaboration has previously reported the first observation of D 0 → h + h − µ + µ − decays, where h is either a pion or a kaon [19]. Charge-conjugate decays are implied throughout, unless stated otherwise. The measured branching fractions are in agreement with SM predictions [7,8]. Selected angular and CP asymmetries were also measured, with results in agreement with the SM and with CP symmetry [20]. However, a complete angular analysis of a rare charm decay is yet to be performed.
This Letter presents the first measurement of the full set of CP -averaged angular observables and CP asymmetries in D 0 → h + h − µ + µ − decays, together with an updated measurement of the CP asymmetry of the total decay rate, defined as where Γ indicates the total decay rate. The measurement uses proton-proton (pp) collision data corresponding to an integrated luminosity of 9 fb −1 , collected by the LHCb experiment at center-of-mass energies of 7 and 8 TeV (run 1) and of 13 TeV (run 2). The analysis is an extension of that reported in Ref. [20]. It uses approximately three times as many signal decays and results for previously measured observables are superseded. Semileptonic D 0 → h + h − µ + µ − decays are described by five independent kinematic variables: the squared invariant masses of the dimuon and dihadron systems, q 2 ≡ m 2 (µ + µ − ) and p 2 ≡ m 2 (h + h − ), and three decay angles θ µ , θ h , φ, (see Fig. S1 in the Supplemental Material [21]). Here, θ µ is the angle between the µ + direction and the direction opposite to that of the D 0 meson in the dimuon rest frame; θ h is the angle between the h + direction and the direction opposite to that of the D 0 meson in the dihadron rest frame; and φ is the angle between the two planes formed by the dimuon and the dihadron systems in the rest frame of the D 0 meson [8,20,21]. In contrast to Ref.
[20], the same definition of the angles is kept for D 0 and D 0 mesons. Following Ref. [8] and defining Ω ≡ (cos θ µ , cos θ h , φ), the differential decay rate is expressed as the sum of nine angular coefficients I 1−9 that depend on q 2 , p 2 , and cos θ h , multiplied by the terms c 1 = 1, c 2 = cos 2θ µ , c 3 = sin 2 θ µ cos 2φ, c 4 = sin 2θ µ cos φ, c 5 = sin θ µ cos φ, c 6 = cos θ µ , c 7 = sin θ µ sin φ, c 8 = sin 2θ µ sin φ, and c 9 = sin 2 θ µ sin 2φ, as Piecewise integration ranges in φ and cos θ µ can be defined such that the coefficients I 2−9 are expressed as angular asymmetries. For example, the coefficient I 2 is obtained as Corresponding integration ranges to obtain I 3−9 are reported in Refs. [8,21]. The coefficients I 2,3,4,7 are even under CP transformations, while I 5, 6,8,9 are CP odd. Since the term c 1 has no dependence on the decay angles, I 1 provides only a normalization factor and is not considered in this Letter. This analysis measures the normalized observables I 2−9 defined as where Γ is the decay rate in the considered region of dimuon mass. The integration boundaries q 2 min , q 2 max , and p 2 max depend on the dimuon-mass region, where p 2 min = 4m 2 h and m h denotes the hadron mass. The integration in cos θ h is defined to optimize the sensitivity to beyond-SM effects by integrating out contributions from the dominant P-wave resonances in the dihadron system, which further decay into h + h − [8]. Experimentally, the observables are determined by measuring the decay-rate asymmetries of the data split by angular tags defined according to the piecewise integration of the decay rate. As an example, from Eqs. (3) and (4), I 2 is measured as The observables I i , measured separately for D 0 and D 0 mesons, are labeled as I i and I i , respectively. Their CP average S i and asymmetry A i are defined as according to the presence of the known intermediate resonances [20]. For D 0 → π + π − µ + µ − decays the regions are (low-mass) below 525 MeV/c 2 , (η) 525-565 MeV/c 2 , (ρ 0 /ω low) 565-780 MeV/c 2 , (ρ 0 /ω high) 780-950 MeV/c 2 , (φ low) 950-1020 MeV/c 2 , (φ high) 1020-1100 MeV/c 2 , and (high mass) above 1100 MeV/c 2 . For D 0 → K + K − µ + µ − decays, three regions are considered: (low mass) below 525 MeV/c 2 , (η) 525-565 MeV/c 2 , and only one region that combines the low and high ρ 0 /ω region due to limited signal yields. The asymmetries are determined only in m(µ + µ − ) regions where a significant signal yield was previously observed [19]. No measurement is performed in the η region of both channels and in the high-mass region of D 0 → π + π − µ + µ − . The integrated measurement includes candidates from all m(µ + µ − ) intervals. For each m(µ + µ − ) region, the kinematically allowed m(h + h − ) range up to a maximum of 1200 MeV/c 2 , is considered. To avoid potential experimenter bias on the measured quantities, the observables were shifted by an unknown value during the development of the analysis and examined only after the analysis procedure was finalized. The LHCb detector is a single-arm forward spectrometer described in detail in Ref. [22]. It provides high-precision tracking and good particle identification over a large range in momentum [23]. Simulation [24][25][26][27][28][29][30][31] is used to optimize the selection and to estimate variations of the reconstruction and selection efficiency across the decay phase space. Corrections to account for mismodeling of the charged-particle multiplicity of the events and of the particle-identification performance are applied using control channels in the data [32,33].
Events are selected online by a trigger that consists of a hardware stage, based on information from the muon systems, followed by a software stage, based on charged tracks that are displaced from any primary pp-interaction vertex (PV). A subsequent software trigger exploits a full event reconstruction [34] Online selection requirements that have changed over data-taking periods are equalized in the off-line selection, which follows closely that of Ref. [20].
Candidate D 0 mesons are constructed off-line by combining four charged tracks, each having momentum larger than 3000 MeV/c and the momentum component transverse to the beam direction p T > 300 MeV/c, that form a good-quality secondary vertex (SV) significantly displaced from any PV in the event. Two oppositely charged particles are required to be identified as muons and two as either pions or kaons. The D 0 candidates are required to have invariant mass in the range 1820 < m(h + h − µ + µ − ) < 1940 MeV/c 2 and to be consistent with originating from the associated PV, defined as the PV with respect to which the D 0 candidate has the lowest impact-parameter significance. The D 0 momentum is required to be aligned with the vector connecting the PV and the SV. The mass of the dihadron system is required to be less than 1200 MeV/c 2 . The D 0 candidates are combined with low-momentum charged pions having p T > 120 MeV/c, denoted as soft pions in the following, to form D * + candidates. The D * + decay vertex is required to coincide with the position of the associated PV. The difference between the masses of the D * + and D 0 candidates, ∆m, must be within |∆m − 145.4 MeV/c 2 | < 0.6 MeV/c 2 , corresponding to approximately ±2 standard deviations in ∆m resolution around the known value [35].
To suppress the combinatorial background formed with randomly associated tracks that accidentally fulfil the selection requirements, a boosted decision tree (BDT) algorithm [36,37] with gradient boosting [38] is employed. Simulated decays and data candidates from the sideband region m(h + h − µ + µ − ) > 1890 MeV/c 2 are used as signal and background proxies, respectively. To have an unbiased estimate of the BDT performance, a cross validation is performed. The training samples are randomly split into two halves and the BDT classifier is applied to the subsample that has not been used in the training. Separate classifiers are trained for D 0 → π + π − µ + µ − and D 0 → K + K − µ + µ − decays and for run 1 and run 2 data samples to account for differences in decay kinematics and data-taking conditions, respectively. The variables used in the training are momentum and p T of the soft pion, the largest distance of closest approach of the D 0 decay-product trajectories, the angle between the D 0 momentum and the vector connecting the PV and the SV, the fit quality of the SV and its spatial separation from the PV. Purely hadronic decays of the form D 0 → h + h − π + π − with two pions wrongly identified as muons are further reduced by requirements on muon identification [39,40]. The optimal working points of the BDT output selection and muon-identification thresholds are determined simultaneously by maximizing the quantity S/ √ S + B, where S and B are the signal and background yields, respectively, determined from the data in the signal region defined as 1840 < m(h + h − µ + µ − ) < 1890 MeV/c 2 . In the approximately 0.5% of events where multiple D 0 → π + π − µ + µ − candidates are reconstructed after the full selection, only one is kept at random. No multiple-candidate events are found for The m(h + h − µ + µ − ) distributions for selected candidates are shown in Fig. 1. Unbinned maximum-likelihood fits to these distributions yield 3579 ± 71 D 0 → π + π − µ + µ − and 318 ± 19 D 0 → K + K − µ + µ − signal decays. The signal probability density function (PDF) is described by a Hypatia distribution [41] with parameters fixed from simulation, apart from two factors scaling the width and mean of the distribution to account for data-simulation differences. Misidentified hadronic decays are described by a Johnson S U distribution [42] with parameters fixed from a fit to high-yield data samples of D 0 → h + h − π + π − decays with muon-mass hypothesis assigned to two pions and muon-identification criteria applied only to one of them. The combinatorial background is described by an exponential function with shape fixed from a fit to the candidates satisfying |∆m − 145.4 MeV/c 2 | > 2 MeV/c 2 , ∆m < 170 MeV/c 2 , and The LHCb detector geometry, signal reconstruction, and selection requirements result in nonuniform efficiency across the five-dimensional phase space of the decays defined by p 2 , q 2 , θ µ , θ h , and φ. The efficiency variations are corrected using a method developed in Refs. [20,43]. A BDT classifier with gradient boosting [36][37][38] is used to identify differences in the decay kinematics before and after reconstruction and selection. The BDT classifier is trained on simulation using the five-dimensional phase-space variables as input. From the classifier output, per-event candidate weights are derived that correspond to the inverse efficiency. To account for the different detector conditions, separate classifiers are trained for run 1 and run 2 data samples. As a consequence of the weighting, the effective statistical power of the To determine the CP asymmetry A CP , the raw asymmetry in D 0 -and D 0 -signal yields, A raw CP , is corrected for O(1%) nuisance asymmetries: differences in the production cross section of D * + and D * − mesons, A P , and in the detection efficiencies of positively and negatively charged soft pions, A D . The raw asymmetry for decays to the final state f is approximated as A raw [44]. In this procedure, the two-dimensional distribution of momentum and pseudorapidity of D * + candidates in the control samples is equalized to that of the signal decays to account for differences in decay kinematics. Since the angular observables are measured independently for D 0 and D 0 mesons, no correction for nuisance asymmetries is needed.
Each angular observable (or A raw CP ) is determined, independently in each dimuon-mass region, through simultaneous unbinned maximum-likelihood fits to the efficiency-corrected m(h + h − µ + µ − ) distributions of candidates split according to the angular tag (or D 0 -meson flavor). The fits use the same model as described earlier, but with PDFs determined independently in each dimuon-mass region. The same PDFs are assumed when fitting the subsamples split by the tags, except for the measurement of I 2 , where the mass shape of misidentified hadronic decays depends on the angular tag. The yield and angular observable (or A raw CP ) of the three fit components (signal, misidentified, and combinatorial background) are the only floating parameters. The results for S i , A i , and A CP , including both statistical and systematic uncertainties added in quadrature, are reported in Fig. 2. In general, the null-test observables show agreement with the SM predictions. A tabulated version is given in the Supplemental Material [21] together with the correlations between the observables, estimated using a bootstrapping technique [45].
Systematic uncertainties are typically between 10% and 50% of the statistical uncertainty, depending on the observable and the dimuon-mass region. These arise from the following sources: the model used in the mass fits; neglected background from partially reconstructed D + s mesons, from D * + candidates made of D 0 mesons combined with unrelated soft pions, and from D * + candidates originating from decays of b-flavored hadrons; uncertainties in the estimation of the efficiency correction; the accuracy of the correction for nuisance charge asymmetries; the finite resolution of the angular variables. The leading systematic uncertainties are those related to the efficiency correction. These include residual biases on the observables due to efficiency variations that are not fully No measurement is performed in the regions indicated by the vertical gray bands. The horizontal bands correspond to the measurements integrated in the dimuon mass, including candidates from all m(µ + µ − ) ranges. The high-mass region of D 0 → π + π − µ + µ − extends to 1590.5 MeV/c 2 and has been truncated on the plots for a clearer visualization of the other regions.
accounted for by the correction, the uncertainty coming from the limited size of the simulation sample, and effects due to potential residual differences between the data and simulation. The systematic uncertainties due to the efficiency correction are evaluated by repeating the analysis on either fully simulated decays or on simplified simulations that mimic the presence of data-simulation differences. The generation of the simplified simulations reproduces the m(µ + µ − ) and m(h + h − ) distributions observed in the data and is performed exploiting multithreaded architectures using the Hydra library [46,47].
The analysis procedure is validated using the more abundant D 0 → K − π + µ + µ − decay in the dimuon-mass range 675 < m(µ + µ − ) < 875 MeV/c 2 , where the contribution from the ρ 0 /ω → µ + µ − decay is dominant. The decay is not sensitive to FCNC processes and is dominated by the SM tree-level amplitude. The analysis measures the angular observables serving as SM null test ( A 2−9 and S 5−7 ) to be consistent with zero within approximately 1%. As a further cross-check, the analysis is repeated on disjoint subsamples of the data selected according to criteria such as the magnetic-field orientation, which is reversed periodically during data taking; the number of PVs in the event; the transverse momentum of the D * + and soft-pion candidates; and the minimum distance of the D 0 meson to the PV. The resulting variations of the measured observables are as expected according to statistical variations.
In summary, a measurement of the full set of CP -averaged angular observables and their CP asymmetries in D 0 → π + π − µ + µ − and D 0 → K + K − µ + µ − decays is reported, together with an updated measurement of A CP . The analysis uses pp collision data collected with the LHCb detector at center-of-mass energies of 7, 8, and 13 TeV, corresponding to an integrated luminosity of 9 fb −1 . This is the first full angular analysis of a rare charm decay ever performed. The measured null-test observables A CP , S 5−7 and A 2−9 are in agreement with the SM null hypothesis with overall p values of 79% (0.8%) for D 0 → π + π − µ + µ − (D 0 → K + K − µ + µ − ) decays, corresponding to 0.3 (2.7) Gaussian standard deviations. These measurements will help constraining the parameters space of physics models extending the SM.
The variable cos θ µ (cos θ h ) is the cosine of the angle between the momentum of the positive muon (hadron) in the rest frame of the dimuon (dihadron) system with respect to the dimuon (dihadron) flight direction as seen from the rest frame of the D 0 candidate (see Fig. S1): cos θ µ = e µµ · e µ + , cos θ h = e hh · e h + .
Here, e kk (k = µ, h) is the unit vector along the momentum of the dimuon (dihadron) system in the rest frame of the D 0 meson and e k + is the unit vector along the momentum of the positively charged muon (hadron) in the dimuon (dihadron) rest frame. The angle φ is the angle between the two decay planes of the dimuon and dihadron systems, defined by: where n kk = e k + × e k − is defined as the unit vector perpendicular to the decay plane spanned by the two muons (or the two hadrons). The angles are defined in −1 ≤ cos θ h ≤ 1, −1 ≤ cos θ µ ≤ 1 and −π ≤ φ ≤ π. In contrast to Ref.

(S8)
The observables I i , measured separately for D 0 and D 0 mesons, are labelled as I i and I i , respectively. The observables reported in the Letter are the CP averages, S i , and asymmetries, A i , defined as for the CP -even (CP -odd) coefficients I 2,3,4,7 ( I 5, 6,8,9 ). The full set of measured CP -averaged angular observables for both decay modes is reported in Table S2; their CP asymmetries are reported in Table S3; the CP asymmetry in the decay rates A CP , are reported in Table S1.
−2.5 ± 6.8 ± 0.6 Full range −2.3 ± 6.3 ± 0.6 reported in the first column. The first uncertainty is statistical and the second systematic.