Measurement of $CP$ asymmetry in $D^0 \to K^0_S K^0_S$ decays

A measurement of the $CP$ asymmetry in $D^0 \to K^0_S K^0_S$ decays is reported, based on a data sample of proton-proton collisions collected by the LHCb experiment from 2015 to 2018, corresponding to an integrated luminosity of 6 fb$^{-1}$. The flavor of the $D^0$ candidate is determined using the charge of the $D^{*\pm}$ meson, from which the decay is required to originate. The $D^0 \to K^+ K^-$ decay is used as a calibration channel. The time-integrated $CP$ asymmetry for the $D^0 \to K^0_S K^0_S$ mode is measured to be: $$A^{CP}(D^0 \to K^0_S K^0_S) = (-3.1\pm 1.2\pm 0.4 \pm 0.2)\%, $$ where the first uncertainty is statistical, the second is systematic, and the third is due to the uncertainty on the $CP$ asymmetry of the calibration channel. This is the most precise determination of this quantity to date.

The existence of charge-parity violation (CP V ) effects in charm hadrons has recently been established [1], which constitutes the only evidence of CP V in the up-quark (u, c or t) sector. However, current experimental evidence is limited to a single observable, ∆A CP , the difference between the CP asymmetry in D 0 → K + K − and D 0 → π + π − decays. Given the uncertainties in the theoretical predictions on CP V in the charm-quark sector within the Standard Model (SM), it is currently not possible to reach a definitive conclusion about their ability to explain the data [2,3]. Further measurements in charm-hadron decays are crucial to shed light on CP V phenomenology. This could involve dynamics beyond the SM, which is not constrained to be the same as in the down-quark (d, s or b) sector and could enter the amplitudes via loop contributions, affecting observables in a detectable way.
Among many charm-hadron decay modes in which CP V could manifest, the D 0 → K 0 S K 0 S mode is a promising target because of the expected size of the effect. Its CP asymmetry is defined by where Γ is the decay width of the D 0 (D 0 ) meson, and it can be larger than in other channels, up to the percent level in the SM [4][5][6][7][8]. In fact, only amplitudes proceeding via loop-suppressed and tree-level exchange diagrams, which vanish in the flavor-SU(3) limit, contribute to this decay, and they are similar in size. Their interference could therefore result in a detectable CP asymmetry. In addition, the CP asymmetry in the D 0 → K 0 S K 0 S decay is sensitive to a different mix of amplitudes compared to D 0 → K + K − and D 0 → π + π − decays. Therefore, measuring A CP (K 0 S K 0 S ) provides independent information which can help to elucidate the mechanisms of CP V in charm hadron decays.
The current world average for the time-integrated CP asymmetry is A CP (K 0 S K 0 S ) = (0.4 ± 1.4)% [9], the precision of which is still insufficient for detecting an effect. In this work, a new measurement of this quantity, performed with proton-proton (pp) collisions collected from 2015 to 2018 (Run 2) by the LHCb experiment at the LHC at CERN, is reported. Data collected at a center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 6 fb −1 , are used. The subsample of data collected during 2015 and 2016 has been analyzed previously [10]. In the present work a number of improvements in the analysis leads to a more efficient selection and a sizeable improvement in sensitivity.
The measurement of A CP (K 0 S K 0 S ) requires knowledge of the D 0 flavor at production. A sample of flavor-tagged D 0 → K 0 S K 0 S decays is obtained by selecting only D 0 mesons that originate from D * + → D 0 π + decays. 1 The charge of the pion (tagging pion) in this decay identifies the flavor of the accompanying D 0 meson. While D 0 oscillations may cause some of them to change flavor before decaying, this is a small effect in comparison to the resolution of the current measurement and is not considered further.
The decay widths Γ in Eq. 1 are related to the number of observed candidates N by where σ(D * ± ) are the production cross sections and ε ± the detection probabilities for D * ± decays. Both factors are charge asymmetric, due to the D * ± production asymmetry arising from the hadronization of charm quarks in pp collisions, and to asymmetries in the geometry and response of the detector. Because of the charge-symmetric final state of the D 0 decay, the only source of detection asymmetry comes from the tagging pion. All these quantities are calibrated using a large sample of D 0 → K + K − decays, for which A CP is known with much greater precision than in the D 0 → K 0 S K 0 S decay mode [11]. The LHCb detector is a single-arm forward spectrometer designed for the study of particles containing b or c quarks, as described in detail in Refs. [12,13]. The LHCb tracking system exploits a dipole magnet to measure the momentum of charged particles, and it consists of a silicon-strip vertex detector surrounding the pp interaction region [14] and tracking stations placed upstream and downstream of the magnet. The magneticfield polarity is reversed periodically during data taking, alternating between pointing upwards (MagUp) and downwards (MagDown), to mitigate the differences of reconstruction efficiencies of particles with opposite charges. The online event selection is performed by a trigger [15], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction and selection in two steps. The simulated data used in this analysis are produced using the software described in Refs. [16][17][18][19][20].
For an event to be considered in this analysis, the hardware-trigger decision is required to have been based either on the transverse energy deposited in the hadronic calorimeter by a particle associated with the D 0 candidate, or on signatures not associated with the D * + candidate, such as a muon with high transverse momentum, or a large transverse-energy deposit in the electromagnetic or hadronic calorimeter from the rest of the event. Similarly, the first software stage decision is required to have been based either on the presence of a single track associated with the D 0 decay, with sufficient transverse momentum and impact parameter relative to the primary pp collision vertex (PV), or on track-related signatures independent of the D * + candidate. This ensures that the trigger decision is independent of the tagging pion, preventing a possible charge-dependent bias which would be difficult to measure and correct accurately.
At the second software stage, a complete reconstruction of the decay chain is performed, requiring a pair of K 0 S mesons compatible with the decay of a D 0 particle, and the presence of an additional pion, compatible with originating from a D * + → D 0 π + decay. Candidate K 0 S → π + π − decays are classified in two categories: those in which the decay occurred early enough for the pions to be reconstructed in the vertex detector ("long", indicated in the following by an "L" label); and those decaying later, such that segments of pion tracks can only be formed in the silicon-strip tracker upstream of the magnet, and the tracking stations after the magnet ("downstream", indicated in the following by a "D" label). The geometric acceptance of the LHCb detector for typical K 0 S momenta is about a factor of 2 larger for D-type versus L-type decays, but the mass, momentum and vertex resolution of the latter category are better than those of the former. In addition, downstream tracks are only reconstructed in the second-stage software trigger, causing a lower efficiency for D-type decays, as can be deduced from Table 1.
Candidates are separated into three categories: LL, LD, or DD, according to the types of the two K 0 S mesons associated to the D 0 meson. These are analyzed separately, and the results are combined only at the end. To prevent any experimenters' bias, the flavor of D 0 candidates was not examined until the analysis was finalized.
Some of the D * + mesons are not promptly produced in the primary pp interaction, but rather come from the decay of a beauty-flavored hadron. These secondary decays are affected by different production asymmetries with respect to prompt decays, which is estimated at the 1% level. No attempt is made in this analysis to reject secondary decays, as selection requirements based on pointing observables are less effective in D 0 mesons reconstructed with displaced vertices, such as D 0 → K 0 S K 0 S , compared to other charm hadron decays. The spurious asymmetry introduced by this contribution is instead canceled by use of a dedicated calibration sample. For this purpose, a sample of D 0 → K + K − decays produced in D * + decays is selected in a way to ensure it contains prompt and secondary decays in the same proportions as the D 0 → K 0 S K 0 S sample. Its selection, starting from the first trigger level, refrains from any requirements for which the effect on the secondary-to-prompt ratio cannot reliably reproduce that in the signal sample. This requires, among other things, avoiding selections based on observables sensitive to the location of the D * or D 0 vertices, which have different resolutions in the two samples [21].
is the invariant mass of the D * + → D 0 π + candidate. This sample may be contaminated by other physics processes. Amongst them, partially reconstructed D 0 decays coming from D * + mesons, such as D 0 → K 0 S K 0 S π 0 , can also peak in ∆m. Some non-D 0 decays can also contaminate the sample, the main contributor being D + s → K 0 S K 0 S π + decays, where the K 0 S K 0 S pair is incorrectly associated with a D 0 decay. These background components are effectively suppressed by selecting only candidates with m(D 0 ) within ±20(30) MeV/c 2 of the known D 0 mass [9] in LL/LD (DD) samples.
A background source which is difficult to eliminate completely is the abundant D 0 →K 0 S π + π − decay, with the π + π − pair mimicking a K 0 S . After a loose selection based on the flight distance of the K 0 S , it is statistically separated from the signal through a simultaneous fit of three observables: the ∆m and the invariant masses, m(K 0 S ), of the two K 0 S candidates. To maximize the sensitivity of the measurement, the LL and LD subsamples are further split according to whether the D * + candidate is compatible or incompatible with having originated from the PV. Compatibility is defined by a fixed threshold on the goodness of fit of the complete decay chain [22], when constrained to be originating from the PV. For PV-compatible candidates, independent of their true nature, a factor of 2 better mass resolution is obtained by constraining the origin vertex to coincide with the PV. Therefore, separate analyses of the two subsamples, followed by the combination of the results, provides optimal exploitation of the available data. This is not done for the DD sample due to its limited signal yield. While this requirement might in principle have slightly different effects on the secondary fraction of signal and calibration sample, differences are negligible, as it has been explicitly verified on data [21].
The data are further split according to a classifier sensitive to the signal-to-background ratio. The classifier combines a number of track-related observables, including track and vertex quality, transverse momenta of K 0 S and D 0 candidates, helicity angles of the K 0 S and D 0 decays, and particle identification parameters of the D 0 final state particles. The kinematics of the tagging pion is excluded to avoid introducing possible charge asymmetry biases. All these observables are combined in a k-nearest-neighbors (kNN) classifier [23], trained using a simulated signal sample, and data from the D 0 mass sidebands for background. The resulting discriminant is used to split each sample in three categories (low, medium and high purity); the low-purity class contains very little signal and is removed. The thresholds separating the three categories are numerically optimized in a  way to achieve the best combined resolution on A CP [21].
Finally, 2015-2016 data are analyzed separately from 2017-2018 data, due to different trigger conditions. In addition, DD candidates were not collected in 2015-2016 data taking. In about 10% of cases, multiple D * + candidates are found in the same event. This is mainly due to D 0 candidates that are associated with multiple tagging pions. In these cases, one D * + candidate is randomly selected for the analysis.
Distributions of the ∆m variable for some representative subsamples are shown in Fig. 1. An unbinned maximum-likelihood fit is performed to the joint distribution of ∆m and the two m(K 0 S ) observables, simultaneously to candidates of both flavors to obtain the value of A CP . The total probability density function is parameterized by the sum of eight components: the signal component, peaking in the three observables, and seven additional components, each describing a specific background source. This includes D 0 → K 0 S π + π − decays, which peak in ∆m and in one m(K 0 S ) distribution but not in the other, and all possible combinations of unrelated particles. The peaking component in the ∆m distribution is described by a Gaussian function in PV-incompatible samples, or a Johnson S U distribution [24] in PV-compatible samples. The peaking component in the m(K 0 S ) distribution is described by the sum of two Gaussian functions with different widths and a common mean, for both L-and D-category K 0 S candidates. The nonpeaking component in the ∆m distribution is described by an empirical threshold function, while in the m(K 0 S ) distribution it is described by a first-and second-order Chebyshev polynomial for L-and D-category K 0 S candidates, respectively. In each subsample, the parameters defining the signal and background probability density functions are shared between flavors, while the normalization of each component is allowed to differ.
Each candidate participating in the fit is appropriately weighted with the aim of correcting for all spurious asymmetries, with the help of the calibration D 0 → K + K − sample, selected in a way to contain the same proportions of primary and secondary decays. The calibration sample is similarly split between PV-compatible and PV-incompatible categories, and is required to have m(K + K − ) within ±20 MeV/c 2 around the known D 0 mass. In addition, only candidates in a ±1.5 MeV/c 2 range around the ∆m peak are used.
To compute weights, each D 0 → K + K − subsample is classified in categories of observed charge asymmetry, by a kNN classifier based on the space of kinematic parameters of the D 0 candidate. The detector and production asymmetries observed in the calibration sample are independent of the D 0 decay mode, being charge symmetric, and can be used to correct the signal sample for the same effects. This is achieved by weighting each signal candidate in the global fit by a charge-dependent factor where p 0 is the D 0 3-momentum, and n ± C ( p 0 ) is the density of calibration D * ± decays in the p 0 space [21].
To account for a possible dependence of the detection asymmetry on magnet polarity, weights are separately calculated for MagUp and MagDown configurations. Their distributions are shown in Fig. 2, where their difference is clearly visible. To avoid weights affected by large uncertainties, a negligible fraction of candidates having very large weights (greater than 10) are dropped at this stage. The weighting procedure described above is a better alternative to the procedure of binning the D 0 kinematic space and weighting the two samples to have the same distribution. The advantage is a reduction of dimensionality to the single variable of actual relevance, that is the charge asymmetry to be corrected. This reduces the loss in statistical power that occurs when weighting bins of very different populations in the two samples under comparison. In addition, it holds exactly even when the involved asymmetries are not small. This is helpful in the case of the LHCb detector, where some kinematic regions are characterized by sizeable detection asymmetries for charged pions [10]. This affects the tagging pion in 30% of our events.
This weighting method has been extensively checked with both data and simulation. One of the checks is to apply the weighting procedure to half of the D 0 → K + K − sample, using the second half to calculate weights. This has the expected effect of canceling the asymmetry of the sample, that was initially highly significant, (1.3 ± 0.1)%, due to a combination of physics and detector effects.
The systematic uncertainty associated to the limited size of the calibration mode is determined by a bootstrap sampling of the data, and is found to be negligible. The uncertainty due to the presence of residual background in the D 0 → K + K − sample has also been evaluated and found to be a minor effect. Another important source of systematic uncertainty is due to the limited knowledge of the shape of the mass distributions, as the fit function used is purely empirical. Any potential bias in the procedure is evaluated by fitting the model to simulated samples of pseudoexperiments from alternative models. The number and type of peaking components and the behavior of the background at the kinematic threshold are varied, and the largest observed variations are taken as the systematic uncertainty, ranging between 0.4% and 5.6%, depending on the subsample.
A systematic uncertainty is also assessed for possible residual differences in secondary decay fractions between signal and calibration samples. The largest discrepancy between the two samples is the presence of slightly different trigger requirements on the proper decay time of the D 0 candidate. This has been corrected for by increasing the weight of candidates close to the threshold, to emulate the effect of those that have been lost. An uncertainty of this correction is conservatively assessed, that is between 0.1% and 0.2%. Finally, an uncertainty of 0.15% in the input value of A CP (K + K − ) is separately taken into account, measured by the LHCb experiment as A CP (K + K − ) = (0.04 ± 0.12 ± 0.10)% [11].
The results obtained in each subsample are summarized in Table 1. The subsamples corresponding to different kNN classifier ranges are fitted simultaneously with common parameters, so they do not produce separate results. As a check of goodness of fit, a χ 2 has been calculated for each one-dimensional projection of the fit (Fig. 1). All p-values are found to be greater than 0.2.
All partial results in Table 1 are statistically compatible with each other. The weighted average of all measurements using 2015-2016 data, is 2)%, where the first uncertainty is statistical; the second is systematic, obtained by taking the individual contributions as uncorrelated; and the third is from the uncertainty on A CP (K + K − ). This result is compatible with the previously published value based on the same data sample [10], but has a better precision by about 30%, corresponding to an effective doubling of the yields. The sensitivity increase is due to the combined effect of several improvements in the analysis, most notably the new weighting technique, the inclusion of secondary decays, an appropriate categorization of the sample, and the multidimensional likelihood fit.
The asymmetry for the 2017-2018 data is measured to be Treating the systematic uncertainties as uncorrelated between the data taking periods, except for the shape of the fit functions that is considered to be fully correlated, the asymmetry combining the results from the full Run 2 data sample is obtained as 2)%. This measurement supersedes the previous LHCb result [10] and is in agreement with all previous determinations [25][26][27]. It is the most precise measurement of this quantity to date, and it is compatible with no CP asymmetry at the level of 2.4 standard deviations.