Search for time-dependent CP violation in D0 -> K+ K- and D0 -> pi+ pi- decays

A search for time-dependent violation of the charge-parity symmetry in D0→ K+K− and D0→ π+π− decays is performed at the LHCb experiment using proton–proton collision data recorded from 2015 to 2018 at a centre-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 6 fb−1. The D0 meson is required to originate from a D∗(2010)+→ D0π+ decay, such that its flavour at production is identified by the charge of the accompanying pion. The slope of the time-dependent asymmetry of the decay rates of D0 and D0 mesons into the final states under consideration is measured to be ∆YK+K− = (−2.3± 1.5± 0.3)× 10−4, ∆Yπ+π− = (−4.0± 2.8± 0.4)× 10−4, where the first uncertainties are statistical and the second are systematic. These results are compatible with the conservation of the charge-parity symmetry at the level of 2 standard deviations and improve the precision by nearly a factor of two. Published in Phys. Rev. D104 (2021) 072010 © 2021 CERN for the benefit of the LHCb collaboration. CC BY 4.0 licence. †Authors are listed at the end of this paper. ar X iv :2 10 5. 09 88 9v 2 [ he pex ] 2 8 O ct 2 02 1


Introduction
The breaking of the invariance of fundamental interactions under the combined charge conjugation (C) and parity (P ) transformation, commonly named CP violation, is a necessary condition to explain the much larger abundance of matter with respect to antimatter in the universe [1]. Within the standard model (SM) of particle physics, the weak interaction provides a source of CP violation through a single complex phase in the Cabibbo-Kobayashi-Maskawa (CKM) matrix that governs the interaction of quarks with the W boson [2,3]. This CKM paradigm has been tested successfully in the decays of down-type quarks (s or b) in K and B mesons. However, the measured size of CP violation is too small to explain the aforementioned matter-antimatter asymmetry [4], suggesting the existence of additional sources of CP violation beyond the SM.
Hadrons containing charm quarks are the only ones where CP violation and flavourchanging neutral currents (FCNC) involving up-type quarks (u, c or t) can be studied, and provide a unique opportunity to detect new interactions beyond the SM that leave down-type quarks unaffected [5]. Within the SM both CP violation and FCNC for charm hadrons are predicted to be smaller than for kaons and beauty hadrons. The Glashow-Iliopoulos-Maiani mechanism is more effective owing to the smaller mass of the beauty with respect to the top quark and to the smallness of the CKM matrix elements connecting the first two generations of quarks with the third. Furthermore, the contributions from the strange and down quarks cancel in the U -spin limit, where U -spin is the SU (2) subgroup of SU (3) F relating the down and strange quarks. In particular, the combination of CKM matrix elements responsible for CP violation in charm decays in the SM is Im(V cb V * ub /V cs V * us ) ≈ −6 × 10 −4 , corresponding to CP asymmetries typically of the order of 10 −4 to 10 −3 [5].
The LHCb collaboration reported the first observation of CP violation in the decay of D 0 mesons in 2019 [6]. However, theoretical uncertainties on nonperturbative effects of the strong interaction do not allow a rigorous assessment of its compatibility with the SM [5,[7][8][9][10][11]. This has prompted a renewed interest of the theory community in the field [12][13][14][15][16][17][18][19][20]. Complementary searches for time-dependent CP violation in D 0 decays, which has not been observed so far, have the potential to clarify this picture [21].
Cabibbo-suppressed D 0 → f decays, where the final state f = K + K − or π + π − is common to D 0 and D 0 mesons, provide one of the most sensitive tests of time-dependent CP violation through the measurement of the time-dependent asymmetry between the D 0 and D 0 decay rates, where Γ(D 0 → f, t) indicates the decay rate of an initial D 0 meson decaying into the final state f at time t. The dependence of the asymmetry on decay time is due to the oscillation of D 0 into D 0 mesons. This process is parametrised through the mixing parameters x 12 and y 12 , defined as x 12 ≡ 2|M 12 /Γ| and y 12 ≡ |Γ 12 /Γ| [22], where H ≡ M − i 2 Γ is the effective Hamiltonian governing the time evolution of the D 0 -D 0 system and Γ is the average decay width of the mass eigenstates. Since both mixing parameters are smaller than 1% [23][24][25][26][27][28][29], the asymmetry can be expanded to linear order in the mixing parameters as where a d f is the CP asymmetry in the decay, τ D 0 is the lifetime of the D 0 meson, and the ∆Y f parameter is approximately equal to [21] ∆Y f ≈ −x 12 sin φ M f + y 12 a d f .
Here, φ M f is defined as φ M f ≡ arg M 12 A f /Ā f , where A f (Ā f ) indicates the decay amplitude of a D 0 (D 0 ) meson into the final state f . The parameter ∆Y f is approximately equal to the negative of the parameter A f Γ defined as the asymmetry of the effective decay widths of D 0 and D 0 mesons into the final state f , as detailed in Appendix A.
Within the SM, the value of ∆Y f is predicted to be of the order of 10 −5 or less [21,[30][31][32], even though an enhancement up to the level of 10 −4 by nonperturbative effects of the strong interaction is not excluded [21,31]. At the current level of experimental precision, finalstate dependent contributions to ∆Y f can be safely neglected, as detailed in Appendix A. The measurements of ∆Y K + K − and ∆Y π + π − are thus expected to be consistent with each other and, under this assumption, they are collectively denoted as ∆Y. Under the same approximation, the phase φ M f would be equal to a dispersive mixing phase φ M 2 common to all D 0 decays, ∆Y ≈ −x 12 sin φ M 2 [21]. The phase φ M 2 is equal to the phase of M 12 with respect to its ∆U = 2 dominant contribution, hence the subscript "2", and coincides with the mixing phase φ 12 , defined as φ 12 ≡ arg(M 12 /Γ 12 ), in the superweak approximation [21,22,[33][34][35].
Reducing the uncertainty on ∆Y f is also essential to determine the parameter a d K + K − from the measurements of the time-integrated asymmetry of D 0 → K + K − decays [36][37][38][39][40], which is equal to where t K + K − is the average measured decay time, which depends on the experimental environment.
This article presents a new measurement performed using proton-proton (pp) collision data collected by the LHCb experiment at a centre-of-mass energy of 13 TeV in 2015-2018, corresponding to an integrated luminosity of 6 fb −1 . Unlike in Ref. [44], the D 0 meson is required to originate from strong D * (2010) + → D 0 π + tag decays, such that its initial flavour at production is identified by the charge of the tagging pion, π + tag . The inclusion of charge-conjugate processes is implied throughout, except in the discussion of asymmetries. Hereafter the D * (2010) + meson is referred to as D * + .

Measurement overview
The measured raw asymmetry between the number of D 0 and D 0 decays into the final state f at time t, is equal to up to corrections that are of third order in the asymmetries. Here, A det (π + tag ) is the detection asymmetry due to different reconstruction efficiencies of positively and negatively charged tagging pions and A prod (D * + ) is the production asymmetry of D * ± mesons in pp collisions. The measurement of ∆Y f from the slope of A raw (f, t), cf. Eq. (2), is largely insensitive to time-independent asymmetries such as the detection and production asymmetries, which depend only on the kinematics of the particles. However, the requirements used to select and reconstruct the decays introduce correlations between the kinematic variables and the measured decay time of the D 0 meson. This causes an indirect time dependence of the production and detection asymmetries that needs to be accounted for. These nuisance asymmetries are controlled with a precision better than 0.5 × 10 −4 by an equalisation of the kinematics of D * + and D * − candidates, as described in Sect. 5. A further time dependence of A prod (D * + ) arises if the D * + meson is produced in the decay of a B meson instead of in the pp collision. The production asymmetry of these secondary D * + mesons is different from that of D * + mesons originating from the primary pp collision vertex (PV). In addition, the measurement of the decay time of secondary D 0 mesons, which is performed with respect to the PV, is biased towards larger values. The size of this background is assessed based on the distribution of the D 0 impact parameter with respect to its PV and its contribution to the asymmetry is subtracted as detailed in Sect. 6.
Finally, ∆Y f is determined through a χ 2 fit of a linear function to the time-dependent asymmetry, measured in 21 intervals of decay time in the range of 0.45 to 8 τ D 0 .
The analysis method is developed and validated using a sample of right-sign D 0 → K − π + decays -thus named since the charges of the pions from the D * + and D 0 decays have the same sign -which consist mainly of Cabibbo-favoured decays of unmixed D 0 mesons. This control sample has the same topology and kinematic distributions very similar to those of the signal channels, but its dynamical CP asymmetry is known to be smaller than the current experimental uncertainty and can be neglected, as shown in Appendix B. Therefore, the raw asymmetry between the number of D 0 → K − π + and D 0 → K + π − decays is approximately equal to where the right-hand side differs from that of Eq. (6) since it receives no contribution from dynamical CP asymmetry, but contains an additional detection asymmetry from the K − π + final state, A det (K − π + ). This asymmetry is removed by the kinematic equalisation described in Sect. 5, along with the other nuisance asymmetries. The compatibility of the slope of the time-dependent asymmetry of D 0 → K − π + decays, ∆Y K − π + , with zero is thus a useful cross-check of the analysis method. Finally, the D 0 → K − π + sample is also used to estimate the size of the systematic uncertainties that are not expected to differ among the D 0 decay channels, allowing higher precision to be achieved than what would be possible by using the D 0 → K + K − and D 0 → π + π − samples. To avoid experimenter's bias, the analysis method was developed without examining the values of ∆Y h + h − of the signal channels, which were inspected only after the method had been finalised and the systematic uncertainties had been estimated.

LHCb detector
The LHCb detector [46,47] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector 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. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a PV, the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. The magnetic field deflects oppositely charged particles in opposite directions and this leads to detection asymmetries. Therefore, its polarity is reversed around every two weeks throughout the data taking to reduce the effect. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov (RICH) detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
The online event selection is performed by a trigger, which consists of a hardware stage followed by a two-level software stage, which applies a full event reconstruction. At the hardware-trigger stage, events are required to contain a muon with high p T or a hadron, photon or electron with high transverse energy deposited in the calorimeters. For hadrons, the transverse energy threshold is approximately 3.7 GeV. In between the two software stages, an alignment and calibration of the detector is performed in near real-time [48] and updated constants are made available for the trigger, ensuring high-quality tracking and particle identification (PID) information. The excellent performance of the online reconstruction offers the opportunity to perform physics analyses directly using candidates reconstructed at the trigger level [49,50], which the present analysis exploits. The storage of only the triggered candidates enables a reduction in the event size by an order of magnitude.
Simulation is used to estimate the size of the background of secondary D * + mesons from B decays in Sect. 6, and of three-body decays of charm mesons in Sect. 7. In the simulation, pp collisions are generated using Pythia [51] with a specific LHCb configuration [52]. Decays of unstable particles are described by EvtGen [53], in which final-state radiation is generated using Photos [54]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [55] as described in Ref. [56].

Candidate selection
The D * + → D 0 π + tag decay, where the D 0 meson subsequently decays into one of the following h + h − combinations, K − π + , K + K − , or π + π − , is reconstructed at the trigger level. No requirements on the type of hardware-trigger decision are applied, while at least one or both of the tracks from the D 0 decay are required to satisfy the single-or two-track selections of the first-stage software trigger. The former requires the presence of at least one track with high p T and large χ 2 IP with respect to all PVs, where the χ 2 IP is defined as the difference in the vertex-fit χ 2 of a given PV reconstructed with and without the particle being considered. The single-track requirement in the χ 2 IP vs. p T plane changed during the data taking, as illustrated in Fig. 1. In particular, the selection was tighter during 2016. On the other hand, the two-track selection requires the presence of two high-p T tracks forming a good-quality vertex that is significantly displaced from its associated PV, defined as the PV to which the IP of the two-track combination is the smallest. In this case the selection is based on a bonsai boosted decision tree [57] that takes as inputs the χ 2 of the two-track vertex fit, the number of tracks with χ 2 IP > 16, the sum of the p T of the two tracks and the significance of their flight-distance with respect to the associated PV, χ 2 FD . This is defined as the difference in the vertex-fit χ 2 of the PV reconstructed including the two tracks, and the sum of the vertex-fit χ 2 of the PV reconstructed without including them and of the two-track vertex-fit χ 2 . Also for the two-track selection, the requirements employed in 2016, in particular during the data taking with the magnetic field pointing upwards, were tighter with respect to the other years.
The second-stage software trigger combines pairs of oppositely charged tracks with distance of closest approach less than 0.1 mm to form D 0 candidates. Both tracks are required to be of high quality based on the χ 2 per degree of freedom of their track fit (χ 2 /ndf < 3) and on the output of a multivariate classifier trained to identify fake tracks, by combining information from all of the tracking systems. Furthermore, both tracks are required to have p > 5 GeV/c and to have a χ 2 IP with respect to all PVs in the event greater than 4. The tracks are given a pion-or kaon-mass assignment, based on the information from the RICH detectors. The D 0 decay vertex is required to be significantly displaced from the PV, and the angle between the D 0 momentum and the vector connecting the PV and the D 0 decay vertex is required to be less than 1 • . Finally, all remaining good-quality tracks of the event, as described above, which satisfy p > 1 GeV/c and p T > 200 MeV/c, are assigned a pion-mass hypothesis and are combined with the D 0 candidate to form a D * + candidate, the vertex fit of which is required to be of good quality.
In the offline selection, the pseudorapidity of all h + , h − and π + tag tracks is required to lie in the range 2 to 4.2 to exclude candidates that traversed detector material corresponding to more than 0.3 interaction lengths between the pp interaction point and the end of the tracking system, as these candidates exhibit larger detection asymmetries [58]. The D 0 flight distance in the plane transverse to the beam is required to be less than 4 mm to remove D * + candidates produced by hadronic interactions with the detector material, and the z coordinate of the D 0 decay vertex is required to lie within 20 cm from the pp interaction point. 1 The h + h − invariant mass, m(h + h − ), is required to lie in the range [1847.8, 1882.6], [1850.6, 1879.9] and [1846.2, 1884.2] MeV/c 2 for the D 0 → K − π + , D 0 → K + K − and D 0 → π + π − candidates, respectively, corresponding to ±2 times the mass resolution around the known D 0 mass [59]. Finally, to suppress the background from D 0 → K − e + ν e decays, the kaons from the D 0 → K + K − decay are required not to be identified as electrons or positrons, based on the output of a multivariate classifier combining information from all of the detectors. This requirement is applied to both particles to avoid introducing different efficiencies for D 0 and D 0 decays owing to possible PID asymmetries.
In order to improve the resolution on the D 0 decay time, a kinematic fit is performed in which the D * + candidate is required to originate from the associated PV [60]. The resulting average decay-time resolution is 0.11 τ D 0 . At the same time, the resolution on the invariant mass of the D * + candidates, m(D 0 π + tag ), is improved by a factor of two. However, the decay time of D 0 mesons coming from secondary D * + mesons produced in the decay of B mesons is overestimated. The IP of these background D 0 mesons is, in general, greater than zero, contrary to signal candidates, whose IP is equal to zero within the experimental resolution. The background from B-meson decays is suppressed to the 4% level by requiring that the D 0 IP is less than 60 µm and that its decay time is less than 8 τ D 0 . Finally, the D 0 decay time is required to be greater than 0.45 τ D 0 to exclude candidates with low reconstruction efficiency.
After these requirements, around 2.5%, 4.7% and 4.9% of the D 0 → K − π + , D 0 → K + K − and D 0 → π + π − candidates are combined with more than one π + tag candidate to form a D * + candidate. In this case, one D * + candidate per event is selected at random. The distributions of m(D 0 π + tag ) of selected candidates are displayed for the three decay channels in Fig. 2. This quantity is calculated using the known D 0 mass in the determination of the D 0 energy. This choice minimises the impact of the resolution on the invariant mass of the D 0 candidate. The m(D 0 π + tag ) signal window is defined as [2009.2, 2011.3] MeV/c 2 and retains about 96.9% of the signal. The purity within this window is 97.7%, 95.5% and 94.1% for the D 0 → K − π + , D 0 → K + K − and D 0 → π + π − samples, respectively. The residual background is dominated by real D 0 mesons associated with uncorrelated particles and is subtracted by using background candidates in the lateral mass window [2015,2018] MeV/c 2 , weighted with a suitable negative coefficient. The coefficient is determined based on a binned maximum-likelihood fit to the m(D 0 π + tag ) distribution, which relies on an empirical model. In particular, the signal probability density function (PDF) is described by the sum of two Gaussian functions and a Johnson The signal window and the lateral window employed to remove the combinatorial background (grey filled area) are delimited by the vertical dashed lines. Fit projections are overlaid.
where the µ and σ parameters are approximately equal to the mean and standard deviation of the Gaussian-like core, and the δ and γ parameters describe its asymmetric tails. The background PDF, instead, is modelled by the function where m 0 is defined as the sum of the D 0 and π + masses and the small parameters α and β quantify the deviations from a square-root function. The background subtraction is performed without distinguishing between D * + and D * − candidates, but separately in The signal window is delimited by the vertical dashed lines. each decay-time interval. The 21 intervals of decay time, which span the range [0. 45,8] τ D 0 , are chosen to be equally populated, except the last two intervals, which contain half the number of candidates with respect to the others.
The m(h + h − ) distributions after the removal of the m(D 0 π + tag ) background are displayed in Fig. 3. The number of candidates in the signal region is 519, 58 and 18 millions for the K − π + , K + K − and π + π − decay channels, respectively. The number of candidates per integrated luminosity is by a factor of 3.4 larger than that of the measurement with the data collected in 2011-2012 [43], owing to the increased charm-quark production cross-section at the higher centre-of-mass energy [62,63], to the increased trigger rate allowed by the real-time reconstruction of the events [49,50], and to the implementation of the two-track selection in the first stage of the software trigger, which increases the selection efficiency at low decay times.  Candidates/(0.01 rad) The regions of the distributions with asymmetries whose magnitude is larger than 20% are plotted with the same colour as ±20%.
These regions are discarded from the data sample after the kinematic weighting.

Momentum-dependent asymmetries
The data sample is affected by momentum-dependent nuisance asymmetries. The largest of these arise from the π + tag meson and are caused by the vertical magnetic field, which bends oppositely charged particles in opposite directions. For a given magnet polarity, lowmomentum particles of one charge at large or small emission angles in the horizontal plane may be deflected out of the detector or into the uninstrumented LHC beam pipe, whereas particles with the opposite charge are more likely to remain within the acceptance. This is shown in Fig. 4, where the momentum of the π + tag meson is parametrised through its emission angles in the bending and vertical planes, θ x(y) ≡ arctan(p x(y) /p z ), and its curvature in the magnetic field, k ≡ 1/ p 2 x + p 2 z . These asymmetries cancel to a large extent in the average between the samples collected with the magnet polarities pointing upwards (MagUp) and downwards (MagDown), but after averaging smaller residual asymmetries remain including due to right-left misalignment of detector elements and to the nonzero x coordinate of the collision point, to different beam-beam crossing angles for the MagUp and MagDown polarities, and to variations of the detection efficiency over time. Other momentum-dependent asymmetries, which are independent of the magnet polarity, are the D * + production asymmetry and the track-efficiency asymmetry. The latter is caused by the higher occupancy of the detector part downstream of the magnet towards which the negatively charged particles are bent, owing to the large number of electrons produced in the particles interaction with the detector material. Finally, for the D 0 → K − π + decay channel, the asymmetry caused by the different interaction cross-section of positively and negatively charged kaons and pions with matter (with the latter being much smaller) is independent of the magnet polarity.
Since the Q-value of the D * + decay is small with respect to the pion mass, both the magnitude and the direction of the momenta of the D * + , π + tag and D 0 mesons are highly correlated. As a consequence, all aforementioned asymmetries reflect into momentumdependent asymmetries of the D 0 candidate, with all being of similar size. These asymmetries would not bias the measurement of ∆Y h + h − if they did not depend on the D 0 decay time. However, even if the momentum of the D 0 meson is uncorrelated with its decay time, the selection requirements introduce correlations between their measured values. For example, due to the requirement on the χ 2 FD of the D 0 candidate, low decay-time values are measured only if the D 0 momentum is sufficiently large. The largest correlations concern the D 0 transverse momentum, the normalised distribution of which is plotted for each decay-time interval in Fig. 5 (top). The raw asymmetry of the data collected with the MagUp polarity, of the order of 1%, increases as a function of transverse momentum, and correspondingly decreases as a function of decay time, as shown in Fig. 5 (centre) and (bottom). As a result, the dependence on decay time is not linear. The data collected in 2016 present a much larger slope of the time-dependent asymmetry, even if their momentum-dependent asymmetries are similar to those of data collected in 2017, since the correlations induced by the first-stage software trigger during 2016 are larger. The asymmetry slopes for the data collected with the MagDown polarity are considerably smaller, as a result of smaller observed momentum-dependent asymmetries.
These nuisance asymmetries are removed by equalising the kinematics of π + tag and π − tag candidates and of D 0 and D 0 candidates. This is obtained by weighting their kinematic distributions to their average. The weighting is performed with a binned approach in two steps. The first equalises the (θ x , θ y , k) distributions of π + tag and π − tag candidates to remove the largest acceptance and detector asymmetries, and employs 36 intervals in the range [−0.27, 0.27] rad for θ x , 27 intervals in the range [−0.27, 0.27] rad for θ y , and 40 intervals in the range [0, 0.8] c/GeV for k. For each variable, all intervals have the same width. In addition, intervals with fewer than 40 π + tag or π − tag candidates, or where the asymmetry between the number of π + tag and π − tag candidates is greater than 20% in magnitude, are removed by setting the corresponding weights to zero. This avoids weights whose value would be prone to large statistical fluctuations or very different from unity. The effect of these requirements is very similar to the application of the fiducial requirements used to remove phase-space regions characterised by large detector asymmetries in Ref. [6], but removes fewer candidates from the data sample.
Even after this weighting procedure, residual asymmetries of about 0.5% and dependent on the D 0 momentum and pseudorapidity are observed [58]. These asymmetries are removed by the second step of the weighting, which considers the tridimensional distribution of (p T (D 0 ), η(D 0 ), η(π + tag )). The first two variables have the largest correlation with decay time, while η(π + tag ) is included to avoid that the weighting of the D 0 kinematics spoils that of the π + tag meson. This second weighting employs 32 intervals in the range [2,18] GeV/c for p T (D 0 ), 25 intervals in the range [2, 4.5] for η(D 0 ) and 22 intervals in the range [2, 4.2] for η(π + tag ). All intervals in each variable have the same width and limits on the minimum number of candidates and on the maximum asymmetry per interval, as in the first weighting, are applied.
While the impact on the result of the second step of the weighting is smaller than that of the first, the corresponding size of the shift in the ∆Y h + h − is of the same order as that of the final statistical uncertainty. In particular, the second step is essential to remove the asymmetries of the momentum distribution of the D 0 meson. For the K − π + decay channel, these receive a contribution from the detection asymmetries of the K − π + final state, which are not eliminated by a dedicated weighting. They are instead removed by the second step of the baseline weighting, as are the other asymmetries affecting the D 0 momentum.
Since the detection asymmetries and the correlations induced by the trigger depend on the data-taking conditions, the weighting is performed separately in eight subsamples, divided according to the year and to the magnet polarity. Furthermore, since the asymmetries are different between the K − π + final state and the signal decay channels, the weighting is performed independently for each decay channel. The weighting slightly modifies the combined momentum distribution of D 0 and D 0 candidates and, consequently, also the m(D 0 π + tag ) distribution. Therefore, the fits performed to calculate the coefficients to subtract the m(D 0 π + tag ) background are repeated after each step of the weighting and the coefficients are updated accordingly.
The measured values of ∆Y h + h − for the three decay channels are displayed in Fig. 6 before and after the kinematic weighting. After the weighting, the time dependence of the asymmetry in each sample is well described by a linear function, as confirmed by the fact that all fits have a χ 2 /ndf compatible with unity. In addition, the measurements of ∆Y h + h − are compatible among different years and magnet polarities. On the other hand, the compatibility of ∆Y K − π + with zero should be confirmed only after the subtraction of the contribution to the asymmetry of secondary D * + mesons from B decays, which is described in Sect. 6. The agreement among the measured values of ∆Y K − π + in the eight subsamples and the compatibility of their average with zero after the aforementioned subtraction confirm the effectiveness of the weighting, which removes even the larger detection asymmetries of this decay channel with a precision three times better than that of the signal samples.
Due to the correlation between the measured decay time and momentum of the D 0 meson, a possible time-dependent asymmetry due to a nonzero value of ∆Y h + h − would reflect into momentum-dependent asymmetries and would be partially cancelled by the kinematic weighting. This would cause a dilution of the true value of ∆Y h + h − . The size of this dilution is measured by introducing an artificial value of ∆Y K − π + in the D 0 → K − π + raw sample, obtained by filtering the candidates according to an efficiency that changes linearly with decay time, with opposite slopes for D 0 and D 0 candidates. The kinematic weighting is then applied and the measured value of ∆Y K − π + is compared to the introduced one. This procedure is repeated for different values of ∆Y K − π + , up to values as large as 100 times the statistical uncertainty of the final measurement. The dilution is found to have a linear effect on the measured value of ∆Y K − π + , which is equal to (96.9 ± 0.1)% of the introduced one. As a cross-check, the same study is performed also for the signal channels, obtaining compatible dilution factors, although less precise. In Fig. 6 and in the following, the results of all decay channels are corrected to account for this dilution factor, using the value measured in the K − π + channel.

Removal of B decays
The background from B decays produces a biasing contribution to the asymmetry, A(t), even after the removal of the nuisance asymmetries described in Sect. 5 Figure 6: Results of the fits to the time-dependent asymmetry for each subsample, for (top) K − π + , (bottom left) K + K − and (bottom right) π + π − final states. In the y-axis labels, the data-taking year is abbreviated with the last two digits only and the magnet polarity MagUp (MagDown) is abbreviated as "Up" ("Dw"), while "Avg." denotes the weighted average. The fit χ 2 /ndf are reported on the right of each plot. Solid red squares and empty black dots correspond to data before and after the equalisation of the π + tag and D 0 kinematics, respectively, and prior to subtracting the contribution of secondary decays described in Sect. 6. asymmetry is equal to where A sig (t) and A B (t) are the asymmetries of signal and secondary decays from B mesons, and f B (t) is the fraction of secondary decays among the D 0 candidates at given decay time t. Since the flight distance of D 0 candidates is measured with respect to their associated PV, the decay time of the secondary background from B decays is overestimated and the fraction f B (t) increases as a function of decay time. Moreover, the asymmetries of signal and secondary decays differ, mainly because of the different production asymmetries of D * + and B mesons, and of the different asymmetries of their selection efficiency at the hardware-level trigger. As a consequence, secondary decays will introduce a bias on the measurement of ∆Y h + h − . The fraction f B (t) is determined through a binned maximum-likelihood fit to the IP(D 0 ) vs. t(D 0 ) bidimensional distribution of D 0 → K − π + candidates, for D 0 and D 0 samples combined. In the fit, the selection requirement on the IP(D 0 ) range is loosened from [0, 60] to [0,200] µm to increase the discriminating power between signal and secondary decays and have a better handle on the latter category. For the same reason, the m(D 0 π + tag ) signal window is enlarged to [2007.5, 2011.3] MeV/c 2 . In fact, for secondary decays the PV constraint biases the measured value of the angle between the D 0 and π + tag momentum, and consequently the m(D 0 π + tag ) invariant mass, to lower values. The template PDFs are taken from a simplified simulation of signal and secondary decays only, whereas all other particles produced in the pp collision are discarded to minimise the usage of computing resources. To reduce small discrepancies between data and simulation, the kinematics of the D 0 meson is weighted to match that of data [64]. The weighting coefficients are calculated using data with IP(D 0 ) < 60 µm (> 100 µm) for signal (secondary) decays. In the fit, the two-dimensional template PDFs are determined from simulation, and only the time-integrated fraction of secondary decays is left free to vary. The ratio of the fit projections, which are shown for three decay-time intervals in Fig. 7, to data agrees with unity within 10-20%, with the largest discrepancies being due to the accuracy of the simulation of the trigger requirements at low decay times. The impact of these discrepancies, which affect similarly signal and secondary decays and cancel to good extent in the calculation of the fraction f B (t), is estimated as a systematic uncertainty in Sect. 7. The dependence of f B (t) on decay time is displayed for the baseline IP(D 0 ) and m(D 0 π + tag ) requirements in Fig. 8 (left). The fraction increases with decay time from around 2% to 7%, corresponding to a time-integrated value of around 4%.
The difference in asymmetry of secondary and signal decays, entering Eq. (10), is measured from D 0 → K − π + candidates satisfying IP(D 0 ) > 100 µm, where the fraction of secondary decays is about 95%. By construction, the weighting of Sect. 5 sets to zero the asymmetry of the candidates satisfying IP(D 0 ) < 60 µm, which are signal decays apart from the 4% contamination of secondary D * + decays. Thus, the size of the time-integrated asymmetry of signal decays after the kinematic weighting is negligible with respect to that of secondary decays, which is equal to the asymmetry at IP(D 0 ) > 100 µm up to a dilution of around 5%. This asymmetry, which is shown in Fig. 8 (right), is compatible with being independent of decay time and amounts to (2.2 ± 0.4) × 10 −3 . The constant behaviour of the asymmetry difference, A B (t) − A sig (t), is in agreement with expectations. The nuisance time-dependent asymmetries of Sect. 5 cancel to good extent in the difference between secondary and signal decays even before the kinematic weighting, since the kinematics of the two categories of decays are similar. Moreover, both the difference of the production and of the selection asymmetries, where the latter is mainly due to particles other than the D * + decay products responsible for the hardware-trigger decision, are expected to depend weakly on momenta. Therefore, the asymmetry difference is not expected to depend on decay time before the kinematic weighting to first order. Since the weighting does not modify the asymmetry difference to first order, this assumption holds after the kinematic weighting as well. In particular, it is verified explicitly in data that the fraction f B (t) and the asymmetry difference are not changed to first order by the kinematic weighting. Therefore, the order in which the weighting and the subtraction of secondary decays are performed does not affect the results. Finally, the dependence of the production asymmetry on the decay time due to B 0 mixing is found to be negligible within the experimental uncertainty, as is the contribution from CP violation, see Sect. 7.
The asymmetry of signal decays is calculated in each interval of decay time by subtracting the term f B (t)(A B − A sig ) from the measured asymmetry, A(t), using the fitted values of f B (t) and of the time-independent asymmetry difference above. The results for the D 0 → K − π + control sample are plotted in Fig. 9. The shift of the ∆Y K − π + value with respect to that of Fig. 6, where the contribution from B-meson decays is not subtracted, is approximately equal to −0.26 × 10 −4 .
Since no differences are expected in f B (t) and in the asymmetry difference among different D 0 decay channels, their estimates for the K − π + control channel are employed to correct the signal samples as well, to minimise the statistical fluctuations on the values of the shift. The results are reported in Sect. 8.

Systematic uncertainties
The main systematic uncertainties on ∆Y f are due to the subtraction of the combinatorial background under the D * + mass peak, the asymmetry of the time-dependent shifts of the peak position for D * + and D * − mesons, and uncertainties in the subtraction of the contribution of the background from B-meson decays. Minor contributions are related to [%]  limitations in the removal of the nuisance asymmetries described in Sect. 5, as well as to the background of misidentified D-meson decays under the D 0 mass peak. Whenever they are not expected to depend on the decay channel, the systematic uncertainties are evaluated relying on the K − π + final state to minimise the statistical uncertainty on their estimated value.
The removal of the background under the D * + mass peak relies on the assumption that the kinematics and the asymmetry of the background are the same in the signal and in the lateral window used for the background subtraction. A systematic uncertainty on this assumption is assigned by repeating the measurement of ∆Y K − π + using three alternative windows, namely [2004.5, 2008.5] MeV/c 2 , [2013,2015] MeV/c 2 and [2018, 2020] MeV/c 2 . No systematic trends are spotted, and additional studies of the background properties do not reveal any significant differences among the four lateral windows. Therefore, the root mean square of the deviations, 0.10 × 10 −4 , is employed as a conservative estimate of the systematic uncertainty. For the K + K − (π + π − ) channel, instead, the systematic uncertainty is calculated by scaling this value by the ratio of the signal-to-background ratios in the K − π + and in the K + K − (π + π − ) channels, yielding 0.20 × 10 −4 (0.28 × 10 −4 ). Uncertainties in the determination of the coefficient used for the background subtraction can cause a bias. A systematic uncertainty on this effect is estimated by repeating the measurement fitting the combinatorial-background distribution using the alternative PDFs employed in Refs. [6,43] instead of the baseline background model. The maximum deviations from the baseline result, which amount to 0.01 × 10 −4 , 0.04 × 10 −4 and 0.05 × 10 −4 for the K − π + , K + K − and π + π − decay channels, respectively, are assigned as systematic uncertainties. Finally, the impact of possible differences in the background PDF for π + tag and π − tag mesons is estimated with the K − π + channel by repeating the fits separately for D * + and D * − candidates. The fits are performed separately for different years and magnet polarities. The deviation from the baseline result, 0.07 × 10 −4 , is taken as systematic uncertainty. The corresponding systematic uncertainties for the K + K − and π + π − samples are calculated by scaling this value to account for different signal-to-background ratios as before, yielding 0.14 × 10 −4 and 0.19 × 10 −4 , respectively. All the systematic uncertainties on the subtraction of the combinatorial background are summed in quadrature, giving the final values 0.12 × 10 −4 , 0.24 × 10 −4 and 0.34 × 10 −4 for the K − π + , K + K − and π + π − channels, respectively.
The In each time interval, the size of the shift is estimated by comparing the D * + and D * − signal distributions. The shift, which is displayed in Fig. 10, is compatible with zero at small decay times, and increases up to ±2 keV/c 2 at large decay times. The impact of this variation is estimated by repeating the measurement of ∆Y K − π + using a time-and flavour-dependent m(D 0 π + tag ) signal window, defined in each time interval by shifting the baseline window for D * + (D * − ) candidates by plus half (minus half) the measured shift. The deviation of ∆Y K − π + from its baseline value, 0.14 × 10 −4 , is taken as systematic uncertainty for all decay channels.
The subtraction of the contribution of B-meson decays from the asymmetry, see Eq. (10), relies on the correct determination of both their fraction as a function of decay time and of their asymmetry difference with respect to signal decays. The impact of the finite precision of the asymmetry difference is equal to 0.04 × 10 −4 and is taken as systematic uncertainty. The time dependence of the asymmetry difference owing to the B 0 mixing is estimated in simulation, and the obtained template PDF is added to the constant function in the fit in Fig. 8 (right). However, its normalisation is compatible with zero and its impact on the measurement is estimated to be less than 0.04 × 10 −4 . Since this value is equal to the systematic uncertainty on the precision on the asymmetry difference in the constant hypothesis, no additional systematic uncertainty is assigned to avoid double counting. The uncertainty in the determination of the fraction f B as a function of time receives two separate contributions. The first comes from the finite size of the simulation sample used to produce the template PDFs. The second, larger contribution, is due to possible discrepancies of the two-dimensional t(D 0 ) vs. IP(D 0 ) distribution between simulation and data. These are estimated using the data subsample where the D * + meson forms a good-quality vertex with a µ − , which provides a pure sample of B-meson decays. The measured differences between data and simulation are found to be of the same order as those observed in the results of the fit to the t(D 0 ) vs. IP(D 0 ) distribution, whose projections are displayed in Fig. 7. Taking into account the sum of these two contributions, the absolute uncertainty on the fraction at high and low decay times, which drives the impact of B-meson decays on the measurement, is equal to 1.0% and 0.7%, respectively. The corresponding uncertainty on the subtraction of the contribution of the asymmetry of B-meson decays is 0.07 × 10 −4 . The two systematic uncertainties arising due to the uncertainty of the asymmetry and fraction of B-meson decays are summed in quadrature, yielding a total systematic uncertainty equal to 0.08 × 10 −4 .
Another source of background contributing to the systematic uncertainty is that of multibody decays of D mesons, where one daughter particle is not reconstructed. If one of the final-state particles is misidentified, a wrong mass assignment can compensate for the underestimation of the invariant mass due to the unreconstructed particle. For D 0 mesons produced in the decay of a D * + meson, these background contributions appear as a peak in the m(D 0 π + tag ) distribution, albeit with a larger width with respect to the signal. Therefore, unlike pure h + h − combinatorial background, they are not removed by the subtraction of the m(D 0 π + tag ) combinatorial background. The same applies also to D 0 → h + h − decays, where one of the daughter particles is misidentified. Even if these misidentified decays mostly lie outside of the m(h + h − ) signal region, they need to be taken into account to determine correctly the contribution of the other background decays. The same observation applies to D + s → K + K − π + decays, where the D + s meson is produced at the PV and the pion is assigned as π + tag . Since their reconstructed m(h + h − ) and m(D 0 π + tag ) masses are anticorrelated, they are not removed by the subtraction of the m(D 0 π + tag ) combinatorial background, which instead causes the appearance of a dip in their m(h + h − ) distribution.
All the background components are studied using a simplified simulation [65], where the decays of unstable particles are described by EvtGen [53], final-state radiation (FSR) is generated using Photos [54] and the acceptance and the momentum, vertex and IP resolutions are simulated in a parametric way. The same simulation is used to determine the FSR distribution of the signal decays with better precision than what would be possible by using the smaller simulated sample described in Sect. 3. However, while the FSR distribution of signal decays is fixed to the results of the simplified simulation, the signal mass resolution and the PDF tails due to decays in flight of pions and kaons into muons are fixed to those measured in the simulation described in Sect. 3. For all decays, simulated events are weighted to reproduce the effect of the PID requirements on the particles reconstructed as coming from the D 0 final state. The weights are calculated with a data-driven method by employing large calibration samples [66,67], and are parametrised as a function of momentum and pseudorapidity. The background contamination in the signal region is estimated through template fits to the m(h + h − ) data distribution. Only the signal and background yields and the resolution of the signal component are varied, the latter to correct for O(10%) discrepancies in the resolution between data and simulation, whereas the background template PDFs and the signal FSR and tails due to decays in flight are fixed to the simulation results. The results of the fit are displayed in Fig. 11. The fitted ratio of the relative normalisation of the background components with respect to the signal agrees with expectations within 15%, except for D + s → K + K − π + decays, for which the discrepancy is at the 35% level. While the agreement with data is not perfect, the projections capture all the main features of the m(h + h − ) distributions and allow an estimate of the size of the background contamination under the D 0 mass peak with a precision sufficient to assess a systematic uncertainty.
For the K − π + decay channel, the largest contamination in the signal window is due to D 0 → K − + ν decays, where + stands for e + or µ + , and amounts to (2.5 ± 0.1) × 10 −4 of the D 0 → K − π + yield. The time-dependent asymmetry of this background is estimated in the [1750, 1780] MeV/c 2 sideband, after subtracting the contribution from signal decays. The total estimated bias on ∆Y K − π + is less than 0.01×10 −4 . For the K + K − final state, the largest background fractions are (8.2 ± 0.8) × 10 −4 for D 0 → K − π + π 0 , (3.7 ± 0.2) × 10 −4 for D 0 → K − + ν and (2.3 ± 0.2) × 10 −4 for D 0 → K − π + decays. Their asymmetries are estimated in the [1750, 1780] MeV/c 2 and [1920,1970] MeV/c 2 sidebands for the D 0 → K − π + π 0 and D 0 → K − π + decay channels, respectively. Measuring the asymmetry of D 0 → K − + ν decays is particularly challenging owing to the tiny fraction of these decays in data. Therefore, its size is conservatively assigned to the maximum value of the asymmetries measured for all other background channels in all of the D 0 → h + h − decay channels. The total bias on the ∆Y K + K − value due to the background components in the m(K + K − ) signal window is estimated to be 0.06 × 10 −4 . Finally, for the π + π − decay channel the only relevant background contribution is due to D 0 → π − + ν decays, whose fraction in the signal region is (2.5 ± 0.2) × 10 −4 . Their asymmetry difference with respect to the signal is estimated in the same way as for D 0 → K − + ν decays for the K + K − final state, and provokes a bias on ∆Y π + π − less than 0.03 × 10 −4 .
The kinematic equalisation of the momentum distribution of π + tag and π − tag and of D 0 and D 0 mesons is performed through a binned approach. While the choice of the concerned variables is optimised to remove the kinematic asymmetries, the intervals size has to be kept large enough to avoid large statistical fluctuations. Therefore, detector-induced, time-dependent asymmetries might not be completely removed by the kinematic weighting if they vary considerably within the intervals. The size of the residual asymmetries is estimated in the D 0 → K − π + sample by reducing progressively the size of the intervals until the measured value of ∆Y K − π + does not change within the statistical uncertainty. A systematic uncertainty of 0.05 × 10 −4 is estimated as the difference between the value of ∆Y K − π + measured with the baseline scheme and its asymptotic value. As a cross-check of the effectiveness of the kinematic weighting in removing the nuisance asymmetries,   Figure 11: Distributions of m(h + h − ) of (top) K − π + , (left) K + K − and (right) π + π − final states, with the results superimposed. The vertical dashed lines delimit the signal region. The m(h + h − ) template PDF of the D + s → K + K − π + decay has a negative contribution to the left of the known D 0 mass due to the subtraction of the m(D 0 π + tag ) background.
alternative configurations of the kinematic weighting acting on different variables have been employed, including that described in Ref. [68]. The baseline configuration minimises the residual asymmetries of all kinematic variables of the D 0 and π + tag mesons after the weighting. However, all weighting procedures that remove satisfactorily the asymmetries of the D 0 momentum provide measurements of ∆Y K − π + within 0.13 × 10 −4 from the baseline value.
All systematic uncertainties are summarised in Table 1. The slope of the time-dependent asymmetry of the control sample is measured to be ∆Y K − π + = (−0.4 ± 0.5 ± 0.2) × 10 −4 , where the first uncertainty is statistical and the second is systematic, and is compatible with zero as expected. Note, however, that this measurement was not performed blindly. Additional robustness tests are performed to check that the measured value of ∆Y h + h − does not display unexpected dependencies on various observables, including the selections that are satisfied by the D 0 candidate at the hardware and at the first software stage of the trigger; the momentum, the transverse momentum and the pseudorapidity of the D 0 and π + tag mesons; the D 0 flight distance in the plane transverse to the beam; the position of the PV along the beamline; and the number of PVs in the event. No significant dependencies of ∆Y h + h − on any of these variables are found. The measurement is repeated for the signal channels, assigning a zero weight in the weighting procedure of Sect. 5 only to the candidates in the tridimensional-space intervals for which the corresponding intervals of the K − π + sample have fewer than 40 candidates or an asymmetry greater than 20%. In this way, the choice of the zero weights is made independent of the value of ∆Y h + h − . The stability of the measurement is further checked as a function of the threshold of the minimum number of candidates and of the maximum asymmetry per interval. The results of all these tests are compatible with the baseline one within the statistical uncertainty. Finally, possible biases due to the decay-time resolution, approximately 0.11 τ D 0 , are determined in simulation to be less than 0.01 × 10 −4 , and thus are neglected.

Results
The time-dependent asymmetries of the D 0 → K + K − and D 0 → π + π − channels, after the kinematic weighting and the subtraction of the contribution from B-meson decays, are displayed in Fig. 12. Linear fits are superimposed, and the resulting slopes are where the first uncertainties are statistical and the second are systematic. Assuming that all systematic uncertainties are 100% correlated, except those on the m(h + h − ) background, which are taken to be uncorrelated, the difference of ∆Y f between the two final states is equal to and is consistent with zero. Neglecting final-state dependent contributions to ∆Y f , the two values are combined using the best linear unbiased estimator [69,70]. is consistent with zero within two standard deviations, and both its statistical and systematic uncertainties are improved by more than a factor of two with respect to the previous most precise measurement [43].
These results are consistent with no time-dependent CP violation in D 0 → K + K − and D 0 → π + π − decays, and improve by nearly a factor of two on the precision of the previous world average [45]. In particular, they tighten the bounds on the size of the phase φ M 2 , which parametrises dispersive CP -violating contributions to D 0 mixing, by around 35% [58].
where τ is defined as τ ≡ Γt and the coefficients c ± On the contrary, the definition of A CP (t) in Eq. (1) employed in Refs. [41][42][43][44]68] and in the present article is always dominated by the first-order terms, since the CP -even second-order terms cancel in the difference in the numerator. In particular, the coefficient of the linear expansion of A CP (t) in Eq. (2) is equal to ∆Y f up to a multiplicative factor of 4|A f | 2 |Ā f | 2 /(|A f | 2 + |Ā f | 2 ) 2 , whose difference with unity is approximately equal to (a d f ) 2 /2 10 −6 [6,40]. This coefficient has been denoted as −A f Γ in Refs. [41][42][43][44]68], neglecting the 1% correction due to y f CP in Eq. (19). The final-state dependent contributions to ∆Y f in Eq. (16) can be isolated by defining φ M f ≡ φ M 2 +δφ f , where φ M 2 is the intrinsic CP -violating mixing phase of D 0 mesons, defined as the argument of the dispersive mixing amplitude M 12 with respect to its dominant ∆U = 2 component, and δφ f is the relative weak phase of the subleading amplitude responsible for CP violation in the decay with respect to the dominant decay amplitude [21]. By defining δ f the strong phase analogue of δφ f , and using δφ f = −a d f cot δ f , Eq. (16) can be written as where the first term is universal and the second encloses the final-state dependence. The term y 12 |a d f | is estimated to be less than 0.1 × 10 −4 by using available experimental data [6,45] and the minimal assumption that a d K + K − and a d π + π − have opposite signs, which is motivated by U -spin symmetry arguments as well as by experimental evidence [6,27]. The factor x 12 y 12 cot δ f can enhance the dependence on the final state, even though the phase δ f is expected to be of O(1) due to large rescattering at the charm mass scale. On the other hand, the SM predictions for φ M 2 are of the order 2 mrad or less [21,[30][31][32], even though enhancements up to one order of magnitude due to low-energy nonperturbative strong interactions cannot be excluded [21,31].
An alternative parametrisation of CP violation and mixing is based on the explicit expansion of the mass eigenstates of H in terms of the flavour eigenstates, |D 1,2 ≡ p |D 0 ± q D 0 , with |p| 2 + |q| 2 = 1 (CP T invariance is assumed). The corresponding mixing parameters are defined as x ≡ (m 2 − m 1 )/Γ and y ≡ (Γ 2 − Γ 1 )/(2Γ), where m 1,2 and Γ 1,2 are the masses and decay widths of the mass eigenstates. Adopting the convention that |D 1 (|D 2 ) is the approximately CP -odd (CP -even) eigenstate, the following relations hold, x 12 ≈ |x| and y 12 ≈ y, up to corrections quadratic in the CP violation parameter sin φ 12 [21,22,35]. In this parametrisation, the parameter ∆Y f defined in Eq. (16) is equal to where φ λ f is defined as φ λ f ≡ arg[−(qĀ f )/(pA f )]. Neglecting terms of order higher than one in the CP -violation parameters (|q/p| − 1), φ λ f and a d f , Eq. (21) can be written as Finally, the dependence on the final state can be separated from the universal component by defining φ λ f ≡ φ 2 − δφ f , see Ref. [21], where φ 2 is a final-state independent weak phase dubbed φ by the HFLAV collaboration [45] and δφ f is the same as above, obtaining where the ∆Y f parameter is defined as A global fit of the mixing and time-dependent CP -violation parameters to all of the charm measurements except the present one, with the assumption that a d f is zero, 2 provides |∆Y K − π + | < 0.3 × 10 −4 at 90% confidence level [58]. This number is around 60% of the precision of the present measurement of ∆Y K − π + . 2 The asymmetries a d f and a d f are expected to be negligible in the SM since CF and DCS decays are not sensitive to quantum chromodynamics penguin and chromomagnetic dipole operators. Experimentally,