Observation of the mass difference between neutral charm-meson eigenstates

A measurement of mixing and $C\!P$ violation in neutral charm mesons is performed using data reconstructed in proton-proton collisions collected by the LHCb experiment from 2016 to 2018, corresponding to an integrated luminosity of 5.4$\textrm{fb}^{-1}$. A total of $30.6$ million $D^{0}\rightarrow K^{0}_{\mathrm{S}}\pi^{+}\pi^{-}$ decays are analyzed using a method optimized for the measurement of the mass difference between neutral charm-meson eigenstates. Allowing for $C\!P$ violation in mixing and in the interference between mixing and decay, the mass and decay-width differences are measured to be $x_{C\!P} = \left[3.97\pm0.46\textrm{(stat)}\pm 0.29\textrm{(syst)}\right]\times10^{-3}$ and $y_{C\!P} = \left[4.59\pm1.20\textrm{(stat)}\pm0.85\textrm{(syst)}\right]\times10^{-3}$, respectively. The $C\!P$-violating parameters are measured as $\Delta x = \left[-0.27\pm0.18\textrm{(stat)}\pm0.01\textrm{(syst)}\right]\times10^{-3}$ and $\Delta y = \left[0.20\pm0.36\textrm{(stat)}\pm0.13\textrm{(syst)}\right]\times10^{-3}$. This is the first observation of a nonzero mass difference in the $D^0$ meson system, with a significance exceeding seven standard deviations. The data are consistent with $C\!P$ symmetry, and improve existing constraints on the associated parameters.

: Depiction of the interference of mixing and decay if a D 0 and a D 0 meson decay to a common final state f .
Neutral charm mesons propagating freely can change (oscillate) into their own antiparticles, as the mass eigenstates are linear combinations of the flavor eigenstates. These flavor-changing neutral currents do not occur at tree level in the Standard Model (SM) and allow for hypothetical particles of arbitrarily high mass to contribute significantly to the process. This can affect the mixing of mesons and antimesons such that measurements of these processes can probe physics beyond the SM [1][2][3][4].
The mass eigenstates of charm mesons can be written as |D 1,2 ⟩ ≡ p|D 0 ⟩ ± q|D 0 ⟩, where p and q are complex parameters and, in the limit of charge-parity (CP ) symmetry, |D 1 ⟩ (|D 2 ⟩) is defined as the CP even (odd) eigenstate. Mixing of flavor eigenstates is described by the dimensionless parameters x ≡ (m 1 − m 2 )c 2 /Γ and y ≡ (Γ 1 − Γ 2 )/(2Γ), where m 1 (2) and Γ 1(2) are the mass and decay width of the D 1 (2) state, respectively, and Γ is the average decay width [5]. In D 0 and D 0 decays to a common final state, f , CP violation in mixing manifests itself if |q/p| ̸ = 1 or in the interference between mixing and decay if ϕ f ≡ arg(qĀ f /pA f ) ̸ = 0. Here A f (Ā f ) denotes the amplitude of the decay process D 0 → f (D 0 → f ). In the D 0 → K 0 S π + π − decay studied in this Letter, CP violation in the decay (|A f | 2 ̸ = |Ā f | 2 ) is not considered, as in the SM it is negligible for the doubly Cabibbo-suppressed (DCS) and Cabibbo-favored (CF) amplitudes contributing to this process. With this assumption, the CP -violating phase is independent of the final state, ϕ f ≈ ϕ ≈ arg(q/p) [6,7].
Sensitivity to the mixing and CP -violating parameters is offered by the self-conjugate, multibody D 0 → K 0 S π + π − decay [14-18]. Inclusion of the charge-conjugate process is implied unless stated otherwise. This final state is accessible in both D 0 and D 0 decays and leads to interference between the mixing and decay amplitudes, as demonstrated pictorially in Fig. 1. The dynamics of the decay are expressed as a function of two invariant masses following the Dalitz-plot formalism, in which a three-body decay is parametrized by a pair of two-body invariant masses [19,20]. The squared invariant mass m 2 (K 0 S π ± ) is denoted as m 2 ± for D 0 decays and m 2 ∓ for D 0 decays. A mixture of DCS and CF decay amplitudes results in large variations of the strong phase and, with mixing, causes a decay-time evolution of the density of decays across the phase space. A joint analysis of the Dalitz-plot and decay-time distributions may be used to determine the mixing parameters. Splitting the sample by flavor of the charm meson at production probes for CP -violating effects. Usage of multibody decay modes is typically challenging, as it requires knowledge of the variation of the hadronic parameters and excellent control of efficiencies, resolutions, and background effects.
This Letter reports on a measurement of the mixing and CP violation parameters in D 0 → K 0 S π + π − decays using the "bin-flip" method [21], a model-independent approach which obviates the need for detailed models of the efficiency, resolution, and contributing amplitudes. Mixing and CP violation are parametrized by z CP and ∆z, which are defined by z CP ±∆z ≡ − (q/p) ±1 (y+ix). The results are expressed in terms of the CP -even mixing parameters x CP ≡ − Im(z CP ) and y CP ≡ − Re(z CP ), and of the CP -violating differences ∆x ≡ − Im(∆z) and ∆y ≡ − Re(∆z). Conservation of CP symmetry implies x CP = x, y CP = y, and ∆x = ∆y = 0. The method has already been employed by the LHCb collaboration, yielding the single most precise measurement of x CP and ∆x [18].
In the bin-flip method, data are partitioned into disjoint regions (bins) of the Dalitz plot, which are defined to preserve nearly constant strong-phase differences ∆δ(m 2 − , m 2 + ) between the D 0 and D 0 amplitudes within each bin [22]. Two sets of eight bins are formed symmetrically about the m 2 + = m 2 − bisector, as illustrated in Fig. 2. The region satisfying m 2 + > m 2 − , which includes regions dominated by the CF D 0 → K * (892) − π + decay, is given a positive index +b, while the opposite region, where the relative contribution from decays following an oscillation is enhanced, is given a negative index −b. The data are further split into 13 bins of decay time, chosen such that the bins are approximately equally populated. The squared-mass and decay-time resolutions are typically 0.006 GeV 2 /c 4 and 60 fs, respectively, which are smaller than the bin sizes used. Thus, they are neglected and accounted for in the systematic uncertainties.
For each decay-time interval (j), the ratio of the number of decays in each negative Dalitz-plot bin (−b) to its positive counterpart (+b) is measured. The usage of ratios minimizes the need for precise modeling of the efficiency variation across phase space and decay time. For small mixing parameters and CP -conserving decay amplitudes, the expected ratios for initially produced D 0 (D 0 ) mesons, The parameter r b is the value of R bj at t = 0, while X b is the amplitude-weighted strongphase difference between opposing bins. Finally, ⟨t⟩ j (⟨t 2 ⟩ j ) corresponds to the average (squared) decay time in each positive Dalitz-plot region where the mixed contribution is negligible, in units of the D 0 lifetime τ = ℏ/Γ [5], calculated directly from backgroundsubtracted data. The other parameters are determined from a simultaneous fit of the observed R ± bj ratios, in which external information on c b ≡ Re(X b ) and s b ≡ − Im(X b ) [22, 23] is used as a constraint.
Samples of D 0 → K 0 S π + π − decays are reconstructed from proton-proton (pp) collisions collected by the LHCb experiment from 2016 to 2018, corresponding to an integrated luminosity of 5.4 fb −1 . The strong-interaction decay D * + → D 0 π + is used to identify the flavor of the neutral charm meson at production. Throughout this Letter, D * + indicates the D * (2010) + meson and soft pion indicates the pion from its decay. The LHCb detector [24, 25] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks.
Decays of K 0 S → π + π − are reconstructed in two different categories: the first involving K 0 S mesons that decay early enough for the pions to be reconstructed in all tracking detectors; and the second containing K 0 S mesons that decay later such that track segments of the pions cannot be formed in the vertex detector, which surrounds the pp interaction (primary vertex) region, resulting in a worse momentum resolution. The latter category contains more candidates but has slightly worse mass and decay-time resolution as well as larger efficiency variations.
The online event selection consists of a hardware stage, selecting events based on calorimeter and muon detector information, followed by two software stages. In the first software stage, the pion pair from the D 0 decay is required to satisfy criteria on momenta and final-state charged-particle displacements from any primary vertex for at least one pion (one-track) or both together with a vertex quality requirement (two-track). The second software stage fully reconstructs D * + → D 0 π + , D 0 → K 0 S π + π − candidates using further requirements on particle identification, momenta, and track and vertex quality. Specific ranges of displacement and invariant mass are imposed on the reconstructed D 0 and K 0 S candidates. Due to differing efficiencies, the sample is split into four categories, depending on whether or not the K 0 S meson is reconstructed in the VELO and whether or not they satisfy the one-track requirement.
Offline, a kinematic fit constrains the tracks to form vertices according to the decay topology, the K 0 S candidate mass to the known value [5], and the D * + candidate to a primary vertex [26]. In the reconstruction of the Dalitz-plot coordinates, an additional constraint on the D 0 candidate mass to the known value improves the resolution. Charm mesons originating from the decays of b hadrons are suppressed by requiring that the D 0 and soft pion candidates originate from a primary vertex. Candidates are rejected if two of the reconstructed tracks use the same hits in the vertex detector. About 6% of the candidates are from collision events in which multiple candidates are reconstructed, usually by pairing the same D 0 candidate with different soft pions. When this occurs, one candidate is chosen randomly, and the rest are removed from the sample.
Signal yields are determined by fitting the distribution of the mass difference between the D * + and D 0 candidates, denoted as ∆m. The signal probability density function is empirically described by a combination of a Johnson S U distribution [27] and two Gaussian functions, one of which shares a mean with the Johnson S U . The background is dominated by real D 0 decays incorrectly combined with a charged particle not associated with a D * + decay, and is modeled with a smooth phase-space-like model, is the Heaviside step function, m π is the charged-pion mass [5], and α and c are free parameters. Figure 3 shows the ∆m distribution of the entire sample, from which the fit identifies (30.585 ± 0.011) × 10 6 signal decays. This represents a factor of 15 larger yield compared to the previous measurement.
To determine the yields used to form the ratios R ± bj , separate fits are performed for each set of Dalitz-plot and decay-time bins bj. The signal model assumes the same parameters for each pair of positive and negative Dalitz-plot bins, and fixes some parameters from a fit integrated over decay time. Fits are performed independently for D 0 and D 0 candidates, as well as for each of the four data subsamples. The measured signal yields are then corrected for two effects that do not cancel in the ratio: experimentally induced correlations between the phase space and decay time, and charge-dependent efficiencies (detection asymmetries).
Online requirements on the displacement and momenta of the D 0 decay products introduce efficiency variations that are correlated between the phase-space coordinates and the D 0 decay time. The effect depends predominantly on the invariant mass of two pions from the D 0 decay, with the efficiency to reconstruct the candidates at low values decreasing significantly at low D 0 decay times. This can bias the measured yield ratios and produce mixing-like trends. To remove this bias, an approach that estimates the relative efficiencies using data is developed. The Dalitz plot is divided into small, rectangular-like regions formed symmetrically across the bisector. Note that these include the portions above and below the bisector, unlike the bins shown in Fig. 2. In the limit of CP symmetry, the contribution of mixing to such symmetric regions depends only on y CP and the hadronic parameters of the D 0 decay [21]. As oscillations result in a migration of decays from one side of the Dalitz plot to the other, and the regions are symmetric with respect to the bisector, there is no effect from x CP . Given a set of inputs for y CP and the hadronic parameters, the contribution of mixing to the decay-time distributions of these regions can be accounted for, such that the remaining differences between regions come from the efficiency correlations. Relative efficiency maps that align the decay-time distributions in all these symmetric regions can then be calculated. Per-candidate weights assigned by the efficiency maps are integrated over the data in each bin using the sPlot method [28] with ∆m as the discriminating variable. This provides correction factors for each of the fitted signal yields.
In calculating the efficiency maps, the strong phase variation within a Dalitz-plot bin is approximated as constant, such that it can be described by the external inputs (s b ). As y CP and s b are parameters of the fit, the correction maps and corresponding correction factors are calculated for a range of values. The smallness of mixing results in smooth variations of the correction factors for a given Dalitz-plot bin, which allows for precise interpolation between the calculated points with polynomials. These polynomials are then incorporated into the fit as a correction that depends on y CP and s b . The correction is calculated for each yield ratio, but is averaged over the initial flavor of the candidates. The procedure has been validated with pseudoexperiments, and a systematic uncertainty is assigned due to the approximation that s b is constant within a bin.
Corrections are also applied in order to take into account detection asymmetries. Due to utilizing ratios of yields, the analysis is insensitive to detection asymmetries of the K 0 S , as well as the soft pion used to tag the flavor of the candidate. However, the kinematics of the pions produced in the D 0 decay depend on the Dalitz-plot coordinate and D 0 flavor. This can result in asymmetric efficiency variations for D 0 and D 0 candidates that imitate CP violation. The two-track π + π − asymmetry, A det (π + π − ), is determined by measuring detection asymmetries in control samples of D + s → π + π + π − and D + s → ϕπ + decays, in which the ϕ meson is reconstructed through a K + K − pair. A randomly chosen π + in the D + s → π + π + π − decay is paired with the π − to form a proxy for the π + π − pair of interest. The D + s → ϕπ + sample is used to cancel asymmetries induced from the remaining π + , A det (π + ), and other sources, such as the trigger selection, A trigger (D + s ), and the production of D + s and D − s mesons in pp collisions, A prod (D + s ). For asymmetries of O(1%), the raw asymmetries A meas can be approximated as The difference of the two measured asymmetries gives the detection asymmetry of the π + π − pair. The control samples are weighted to match the kinematics of the pions from the D 0 → K 0 S π + π − sample. This weighting is done separately for each Dalitz-plot bin. The detection asymmetries are of the order of 10 −3 and are used as corrections to the measured yields. They are included as constraints in the fit along with the associated covariance matrix ∆V asym describing uncertainties coming from the limited size of the calibration samples.
The mixing parameters are determined by minimizing a least-squares function where the yields N and their measured uncertainties σ are scaled by factors for the correlation removal, C bj , and detection asymmetry correction, . The different subsamples are fitted simultaneously, separated between D 0 and D 0 flavors denoted as + and −, including all decay-time intervals j and Dalitz-plot bins b. The parameters X b are constrained with a Gaussian penalty term using the values X EXT b and covariance matrix V EXT from a combination of CLEO and BESIII measurements [22,23]. In the fit, the parameters r b are determined independently for each subsample, as they are affected by the sample-specific variation of the efficiency over the Dalitz plot [21]. To avoid experimenter's bias, the values of x CP , y CP , ∆x, and ∆y were not examined until the full procedure had been finalized. Figure 4 shows the yield ratios with fit projections overlaid for each of the eight Dalitz-plot bins. Deviations from constant values are due to mixing. The fit projection when x CP is fixed to zero is also included and shows the inability of a nonzero y CP value to produce the deviations on its own. Also shown are the differences of ratios between D 0 and D 0 decays, where a significant slope would indicate CP violation.
Systematic uncertainties are assessed from ensembles of pseudoexperiments. These use the D 0 → K 0 S π + π − model of Ref.
[29] to describe the amplitude at t = 0, and the decay-time dependence is incorporated for a range of values of the mixing and CP violation parameters. Different sources of systematic uncertainty are included, and the effect on the measured parameters evaluated. The dominant systematic uncertainty on the mixing parameters comes from reconstruction and selection effects, and amounts to 0.20 × 10 −3 (0.76 × 10 −3 ) for x CP (y CP ). This includes neglecting the decay-time and m 2 ± resolutions and efficiencies, as well as the correction to remove the efficiency correlations. The most important effect for y CP is the approximation of the strong phase to be constant within each bin in the procedure to remove correlations. Contamination from b-hadron decays contributes 0.20 × 10 −3 (0.15 × 10 −3 ) to the x CP (y CP ) uncertainty. Potential mismodeling in the signal yield fits contributes 0.36 × 10 −3 to the y CP uncertainty. Time-dependent detection asymmetries are present mainly in bins that give the best sensitivity to ∆y, resulting in a systematic uncertainty of 0.12 × 10 −3 .
The consistency of the results is tested by repeating the analysis in subsets of the data, divided according to magnet polarity, trigger and K 0 S category, data-taking period, D * + meson kinematics, and other categories. The largest variation occurs for the value of x CP as a function of D * + meson pseudorapidity, where the compatibility, considering statistical uncertainties only, amounts to a p-value of 1.5%, depending on the details of the sample split, whereas the overall p-value for all x CP variations observed is above 8%. The observed variations of the observables x CP , y CP , ∆x and ∆y are all consistent with statistical fluctuations. where the first uncertainty is statistical and the second systematic. The statistical uncertainty contains a subleading component due to the limited precision of the external measurements of the strong phases and control samples used for the detection asymmetry. This amounts to approximately (0.23, 0.66, 0.04, and 0.08) × 10 −3 for x CP , y CP , ∆x, and ∆y, respectively. The measurements are statistically limited, though the systematic uncertainty on y CP is comparable to the statistical uncertainty. The results are used to form a likelihood function of x, y, |q/p|, and ϕ using a likelihood-ratio ordering that assumes the observed correlations to be independent of the true parameter values [30]. The best fit point is In summary, a measurement of mixing and CP violation in D 0 → K 0 S π + π − decays has been performed with the bin-flip method, using pp collision data collected by the LHCb experiment and corresponding to an integrated luminosity of 5.4 fb −1 . This resulted in the first observation of a nonzero value of the mass difference x of neutral charm meson mass eigenstates with a significance of more than seven standard deviations, and significantly improves limits on mixing-induced CP violation in the charm sector.    [20] E. Fabri, A study of τ -meson decay, Nuovo Cim. 11 (1954) 479.
[23] BESIII collaboration, M. Ablikim et al., Model-independent determination of the relative strong-phase difference between D 0 andD 0 → K 0 S,L π + π − and its impact on the measurement of the CKM angle γ/ϕ 3  1 Supplemental material Table 1 summarizes the measured values along with their uncertainties and correlations.  Table 2 gives the derived values for x, y, q/p and ϕ together with the 95.5% confidence interval. Table 3 shows a summary of the uncertainties in this analysis. Fig. 5 shows the Dalitz plot of the background-subtracted D 0 → K 0 S π + π − candidates used in the analysis. No efficiency corrections are applied. All samples are combined.