Model-independent determination of the relative strong-phase difference between $D^0$ and $\bar{D}^0\rightarrow K^0_{S,L}\pi^+\pi^-$ and its impact on the measurement of the CKM angle $\gamma/\phi_3$

Crucial inputs for a variety of $CP$-violation studies can be determined through the analysis of pairs of quantum-entangled neutral $D$ mesons, which are produced in the decay of the $\psi(3770)$ resonance. The relative strong-phase parameters between $D^0$ and $\bar{D}^0$ in the decays $D^0\rightarrow K^0_{S,L}\pi^+\pi^-$ are studied using 2.93~${\rm fb}^{-1}$ of $e^+e^-$ annihilation data delivered by the BEPCII collider and collected by the BESIII detector at a center-of-mass energy of 3.773 GeV. Results are presented in regions of the phase space of the decay. These are the most precise measurements to date of the strong-phase parameters in $D \to K_{S,L}^0\pi^+\pi^-$ decays. Using these parameters, the associated uncertainty on the Cabibbo-Kobayashi-Maskawa angle $\gamma/\phi_3$ is expected to be between $0.7^\circ$ and $1.2^\circ$, for an analysis using the decay $B^{\pm}\rightarrow DK^{\pm}$, $D\rightarrow K^0_S\pi^+\pi^-$, where $D$ represents a superposition of $D^0$ and $\bar{D^0}$ states. This is a factor of three smaller than that achievable with previous measurements. Furthermore, these results provide valuable input for charm-mixing studies, other measurements of $CP$ violation, and the measurement of strong-phase parameters for other $D$-decay modes.

Crucial inputs for a variety of CP -violation studies can be determined through the analysis of pairs of quantum-entangled neutral D mesons, which are produced in the decay of the ψ(3770) resonance. The relative strong-phase parameters between D 0 andD 0 in the decays D 0 → K 0 S,L π + π − are studied using 2.93 fb −1 of e + e − annihilation data delivered by the BEPCII collider and collected by the BESIII detector at a centerof-mass energy of 3.773 GeV. Results are presented in regions of the phase space of the decay. These are the most precise measurements to date of the strong-phase parameters in D → K 0 S,L π + π − decays. Using these parameters, the associated uncertainty on the Cabibbo-Kobayashi-Maskawa angle γ/φ3 is expected to be between 0.7 • and 1.2 • , for an analysis using the decay B ± → DK ± , D → K 0 S π + π − , where D represents a superposition of D 0 andD 0 states. This is a factor of three smaller than that achievable with previous measurements. Furthermore, these results provide valuable input for charm-mixing studies, other measurements of CP violation, and the measurement of strong-phase parameters for other D-decay modes.

I. INTRODUCTION
The study of quantum-correlated charm-meson pairs produced at threshold allows unique access to hadronic decay properties that are of great interest across a wide range of physics applications. In particular, determination of the strong-phase parameters provides vital input to measurements of the Cabibbo-Kobayashi-Maskawa (CKM) [1] angle γ (also denoted φ 3 ) and other CP -violating observables. The same parameters are required for studies of D 0D0 mixing and CP violation in charm at experiments above threshold. The angle γ is a parameter of the unitarity triangle (UT), which is a geometrical representation of the CKM matrix in the complex plane. Within the standard model (SM) all measurements of unitarity-triangle parameters should be self-consistent. The parameter γ is of particular interest since it is the only angle of the UT that can easily be extracted in tree-level processes, in which the contribution of non-SM effects is expected to be very small [2]. Therefore, a measurement of γ provides a benchmark of the SM with negligible theoretical uncertainties. A precise measurement of γ is an essential ingredient in testing the SM description of CP violation. A comparison between this, direct, measurement of gamma, and the indirect determination coming from the other constraints of the UT is a sensitive probe for new physics.
One of the most sensitive decay channels for measuring γ is B − → DK − , D → K 0 S π + π − [4] where D represents a superposition of D 0 andD 0 mesons. Throughout this paper, charge conjugation is assumed unless otherwise explicitly noted. The amplitude of the B − decay can be written as (1) Here, m 2 + and m 2 − are the squared invariant masses of the K 0 S π + and K 0 S π − pairs from the D 0 → K 0 S π + π − decay, f D (m 2 + , m 2 − )(fD(m 2 + , m 2 − )) is the amplitude of the D 0 (D 0 ) decay to K 0 S π + π − at (m 2 + , m 2 − ) in the Dalitz plot, r B is the ratio of the suppressed amplitude to the favored amplitude, and δ B is the CP -conserving strong-phase difference between them. If the small second-order effects of charm mixing and CP violation [4][5][6][7][8]  . While the strong-phase difference can be inferred from an amplitude model of the decay D 0 → K 0 S π + π − , such an approach introduces model-dependence in the measurement. This property is undesirable as the systematic uncertainty associated with the model is difficult to estimate reliably, since common approaches to amplitude-model building break the optical theorem [9]. Instead, the strong-phase differences may be measured directly in the decays of quantumcorrelated neutral D meson pairs created in the decay of the ψ(3770) resonance [4,7]. This approach ensures a modelindependent [10][11][12][13][14] measurement of γ where the uncertainty in the strong-phase knowledge can be reliably propagated.
Knowledge of the strong-phase difference in D → K 0 S π + π − has important applications beyond the measurement of the angle γ in B ± → DK ± decays. First, this information can be used in γ measurements based on other B decays [12,15]. Second, it can be exploited to provide a model-independent measurement of the CKM angle β through a time-dependent analysis ofB 0 → Dh 0 where h is a light meson [16] and B 0 → Dπ + π − [17]. Finally, D → K 0 S π + π − is also a powerful decay mode for performing precision measurements of oscillation parameters and CP violation in D 0D0 mixing [18][19][20][21]. Again, knowledge of the strong-phase differences allows these measurements to be executed in a model-independent manner [20,21]. The ability to have model-independent results is critical as these mea-surements become increasingly precise with the large data sets that will be analyzed at LHCb and Belle II, over the coming decade.
The strong-phase differences in D → K 0 S π + π − have been studied by the CLEO collaboration using 0.82 fb −1 of data [22,23]. These measurements are limited by their statistical precision and would contribute major uncertainties to the measurements of γ, and mixing and CP violation in the charm sector, anticipated in the near future. The BESIII detector at the BEPCII collider has the largest data sample collected at the ψ(3770) resonance, corresponding to an integrated luminosity of 2.93 fb −1 . Therefore it is possible to substantially improve the knowledge of the strong-phase differences, which will reduce the associated uncertainty when used in other CP violation measurements.
The observables measured in this analysis are the amplitude-weighted average cosine and sine of the strongphase difference for D → K 0 S π + π − and D → K 0 L π + π − in regions of phase space. The paper is organized as follows. In Sec. II the formalism of how the strong-phase information can be accessed is discussed along with the description of the phase space regions. The BESIII detector and the simulated data are described in Sec. III. The event selection is presented in Sec. IV. Sections V and VI describe the measurement of the strong-phase parameters and their systematic uncertainties. The impact of these results on measurements of γ is assessed in Sec. VII. This paper is accompanied by a letter submitted to Physical Review Letters [24].

A. Division of phase space
The analysis of the data is performed in regions of phase space. Measurements are presented in three schemes which are identical to those used in Ref. [23]. All schemes divide the phase space into eight pairs of bins, symmetrically along the m 2 + = m 2 − line. The bins are indexed with i, running from −8 to 8 excluding zero. The bins have a positive index if their position satisfies m 2 + < m 2 − , and the exchange of coordinates (m 2 + , m 2 − ) ↔ (m 2 − , m 2 + ), changes the sign of the bin. The choice of division of the phase space has an impact on the sensitivity of the CP violation measurements that use this strong-phase information as input. The schemes are irregular in shape and are shown in Fig. 1. Detailed information on the choice of these regions is given in Ref. [23]. The scheme denoted "equal binning" defines regions such that the variation in ∆δ D over each bin is minimized, and is based on a model developed on flavor-tagged data [25] to partition the phase space. In the half of the Dalitz plot m 2 + < m 2 − , the ith bin is defined by the condition 2π(i − 3/2)/8 < ∆δ D (m 2 + , m 2 − ) < 2π(i − 1/2)/8. (3) A more sensitive scheme for the measurement of γ, denoted as "optimal binning", takes into account both the model of the D 0 → K 0 S π + π − decay and the expected distribution of D decays arising from the process B − → DK − when determining the bins. This choice improves the sensitivity of γ measurements compared to the equal binning by approximately 10%. The third binning scheme, denoted the "modified optimal binning" is useful in analyzing samples with low yields [12]. Although these three binning schemes are based on the D 0 → K 0 S π + π − model reported in Ref. [25], this procedure does not introduce model-dependence into the analyses that employ the resulting strong-phase measurements. The determination of CP violation parameters will remain unbiased, but they may have a loss in sensitivity with respect to expectation, due to the differences between the model and the true strong-phase variation.

B. Event yields in quantum-correlated data
The interference between the amplitudes of the D 0 and D 0 decays can be parameterized by two quantities c i and s i , which are the amplitude-weighted averages of cos∆δ D and sin∆δ D over each Dalitz plot bin. They are defined as and where F i is the fraction of events found in the ith bin of the flavor-specific decay D 0 → K 0 S π + π − . The ψ(3770) has a C = −1 quantum number and this is conserved in the strong decay in which two neutral D mesons are produced. Hence the two neutral D mesons have an antisymmetric wave function. This also means that the two D mesons do not decay independently of one another.
For example, if one D meson decays to a CP -even eigenstate, e.g. K + K − , then the other D meson is known to be a CP -odd state. The analysis strategy is to use doubletagged events in which both charm mesons are reconstructed. The yield of events in which one meson is flavor-tagged, e.g. through the decay K − e + ν e , and the other decays to . The details of determining K i through using flavor-specific decays are described in Sec. V B.
Considering a pair of decays where one D meson decays to CP eigenstate, referred to as "the tag", and the other D meson decays to the K 0 S π + π − final state, the decay amplitude of the D → K 0 S π + π − decay is given by where f CP ± refers to the CP eigenvalue of the D → K 0 S π + π − decay. It is possible to generalize this expression FIG. 1. The (left) equal ∆δD, (middle) optimal and (right) modified optimal binnings of the D → K 0 S,L π + π − Dalitz plot from Ref. [23]. The color scale represents the absolute value of the bin number |i|.
to include decays where the tag D meson decays to a selfconjugate final state rather than a CP eigenstate, assuming that the CP -even fraction, F CP , is known. The number of events observed in the ith bin, M i , where the tag D meson decays to a self-conjugate final state is then given by where h CP is a normalization factor. The value of F CP is 1 for CP -even tags and 0 for CP -odd tags. This parameterization is valuable since it allows for final states with very high or very low CP -even fractions to be used to provide sensitivity to the c i parameters. A good example of such a decay is the mode D → π + π − π 0 where the fractional CP -even content is measured to be F πππ 0 CP = 0.973 ± 0.017 [26]. However, from Eq. (4), the sign of ∆δ D is undetermined if only the values of c i are known from the CP -tagged D → K 0 S π + π − decay. Important additional information can be gained to determine the s i parameters by studying the Dalitz plot distributions where both D mesons decay to K 0 S π + π − . The amplitude of the ψ(3770) decay is in this case given by where the use of the † symbol differentiates the Dalitz plot coordinates of the two D → K 0 S π + π − decays. The variable M ij is defined as the event yield observed in the ith bin of the first and the jth bin of the second D → K 0 S π + π − Dalitz plot, and is given by where h corr is a normalization factor. Equation (9) is not sensitive to the sign of s i , however, this ambiguity can be resolved using a weak model assumption.
In order to improve the precision of the c i and s i parameters it is useful to increase the possible tags to include D → . Hence, where the D → K 0 L π + π − is used as the signal decay, and the tag is a self-conjugate final state, the observed event yield M i is given by (10) where K i and c i are associated to the D → K 0 L π + π − decay. The event yield M ij , corresponding to the yield of events where the D → K 0 S π + π − decay is observed in the ith bin and the D → K 0 L π + π − decay is observed in the jth bin, is given by where s i is the amplitude-weighted average sine of the strongphase difference for the D → K 0 L π + π − decay. In Eqs. (7), (9), (10) and (11), the normalization factors h Here S CP is the yield of events in which one charm meson is reconstructed as the CP -tag where no requirement is placed on the decay of the other charm meson, and S FT ( ) refers to the analogous quantity summed over flavor-tagged decays that are used in the determination of K ( ) i . The effective efficiency for detecting the D → K 0 S(L) π + π − decay recoiling against the particular CP -tag under consideration, is defined ST is the detection efficiency for finding the CP -tagged candidate, while DT is the efficiency for simultaneously finding the CP -tagged candidate and the signal decay D → K 0 S(L) π + π − . Furthermore, respectively. Note that, as is discussed in Sec. V B, finite detector resolution results in the migration of reconstructed events between Dalitz plot bins. In order to avoid biases arising from these migration effects it is necessary to modify Eqs. (7), (9), (10) and (11)

III. THE BESIII DETECTOR
BEPCII is a double-ring e + e − collider with a center-ofmass energy ranging from 2 to 5 GeV and a design luminosity of 10 33 cm −2 s −1 at a beam energy of 1.89 GeV. The BESIII detector at BEPCII is a cylindrical detector with a solid-angle coverage of 93% of 4π. The detector consists of a helium-gas based main drift chamber (MDC), a plastic scintillator timeof-flight (TOF) system, a CsI(Tl) electromagnetic calorimeter (EMC), a superconducting solenoid providing a 1.0 T magnetic field and a muon counter. The charged-particle momentum resolution is 0.5% at a transverse momentum of 1 GeV/c, and the specific energy loss (dE/dx) resolution is 6% for the electrons from Bhabha scattering. The photon energy resolution in the EMC is 2.5% in the barrel and 5.0% in the end-caps at energies of 1 GeV. The time resolution of the TOF barrel part is 68 ps, while that of the end-cap part is 110 ps. More details about the design and performance of the detector are given in Ref. [27].
A GEANT4-based [28] simulation package, which includes the geometric description of the detector and the detector response, is used to determine signal detection efficiencies and to estimate potential backgrounds. The production of the ψ(3770), initial-state radiation (ISR) production of the ψ(2S) and J/ψ, and the continuum processes e + e − → τ + τ − and e + e − → qq (q = u, d and s) are simulated with the event generator KKMC [29], with the inclusion of ISR effects up to second-order corrections [30]. The final-state radiation effects are simulated via the PHOTOS package [31]. The known decay modes are generated by EVTGEN [32] with the branching fractions (BFs) set to the world average values from the Particle Data Group [33], while the remaining unknown decay modes are modeled by LUNDCHARM [34]. The generation of simulated signals D 0 → K 0 S π + π − and D 0 → K 0 L π + π − is based on the knowledge of isobar resonance amplitudes from the Dalitz plot analysis of D 0 → K 0 S π + π − . For other multibody decay modes the simulated data are also based on amplitude models, where available, or through an estimate of the expected intermediate resonances participating in the decay.

IV. EVENT SELECTION
In order to measure c i , s i , c i and s i , a range of single-tag (ST) and double-tag (DT) samples of D decays are reconstructed. The ST samples are those where the decay products of only one D meson are reconstructed. The DT samples are those where one D meson decays to the signal mode K 0 S π + π − or K 0 L π + π − and the other D meson decays to one of the tag modes listed in Table I. Tag decay modes fall into the categories of flavor, CP eigenstates or mixed-CP . Flavor tags identify the flavor of the decaying meson through a semileptonic decay or a Cabibbo-favored hadronic decay (contamination from doubly-Cabibbo-suppressed (DCS) decays is discussed later). CP eigenstates and mixed-CP tags identify a decay from an initial state which is a superposition of D 0 and D 0 . The D → π + π − π 0 tag is used for the first time to measure the strong-phase parameters in D → K 0 S,L π + π − decays. It has a relatively high BF and selection efficiency resulting in a large increase to the CP -tagged yields. The use of this tag is possible through the knowledge of F CP for this decay [26]. In this paper the D → π + π − π 0 is referred to as a CP -even eigenstate, although its small CP -odd component is always taken into account, as in Eq. (7).
Due to the hermetic nature of the detector it is possible to use missing energy and momentum constraints to infer the presence of the neutrino in the K + e −ν e final state that does not leave a response in the detector. Similarly, the K 0 L meson, which does not decay within the detector, can be inferred by requiring the missing energy and momentum to be consistent with a K 0 L particle. Tag decay modes such as D → K 0 L ω are not included in the analysis as the systematic uncertainty due to the need to estimate their BFs would be larger than the impact on statistical precision brought from the increased CPtag yields. The principles of missing energy and momentum can also be used to increase the selection efficiency in highly sensitive decay modes by only partially reconstructing the D → K 0 S π + π − candidate. The DT combinations that result in two missing particles are not pursued due to the inability to reliably allocate the missing energy and momentum between two missing particles. The ST yields are only measured in decay modes that are fully reconstructable.
In this paper, we use the following selection criteria to reconstruct the ST and DT samples. The charged tracks are required to be well reconstructed in the MDC detector with the polar angle θ satisfying | cos θ| < 0.93. Their distances of the closest approach to the interaction point (IP) are required to be less than 10 cm along the beam direction and less than 1 cm in the perpendicular plane. For tracks originating from K 0 S , their distances of closest approach to the IP are required to be within 20 cm along the beam direction.
To discriminate pions from kaons, the dE/dx and TOF information are used to obtain particle identification (PID) likelihoods for the pion (L π ) and kaon (L K ) hypotheses. Pion and kaon candidates are selected using L π > L K and L K > L π , respectively. To identify the electron, the information measured by the dE/dx, TOF, and EMC are used to construct likelihoods for electron, pion and kaon hypotheses (L e , L π and L K ). The electron candidate must satisfy L e > 0.001 and L e /(L e + L π + L K ) > 0.8. K 0 S mesons are reconstructed from two oppositely charged tracks with an invariant mass within (0.485, 0.510) GeV/c 2 . A fit is applied to constrain these two charged tracks to a common vertex, and the decay vertex is required to be separated from the interaction point by more than twice the standard deviation (σ) of the measured flight distance (L), i.e., L/σ L > 2, in order to suppress the background from pion pairs that do not originate from a K 0 S meson.
Photon candidates are reconstructed from isolated clusters in the EMC in the regions | cos θ| ≤ 0.80 (barrel) and 0.86 ≤ | cos θ| ≤ 0.92 (end cap). The deposited energy of a neutral cluster is required to be larger than 25 (50) MeV in barrel (end cap) region, and the angle between the photon candidate and the nearest charged track must be larger than 10 • . To suppress electronic noise and energy deposits unrelated to the event, the difference between the EMC time and the event start time is required to be within (0, 700) ns. To reconstruct π 0 (η) candidates, the invariant mass of the accepted photon pair is required to be within (0.110, 0.155)[(0.48, 0.58)] GeV/c 2 . To improve the momentum resolution, a kinematic fit is applied to constrain the γγ invariant mass to the nominal π 0 (η) mass [33], and the χ 2 of the kinematic fit is required to be less than 20. The fitted momenta of the π 0 (η) are used in the further analysis. When reconstructing η candidates decaying through η → π + π − π 0 , it is required that their invariant masses be within (0.530, 0.655) GeV/c 2 . Similarly, ω candidates are selected by requiring the invariant mass of π + π − π 0 to be within (0.750, 0.820) GeV/c 2 . The decay modes η → π + π − η and η → γπ + π − are used to reconstruct η mesons, with the invariant masses of the π + π − η and γπ + π − required to be within (0.940, 0.976) and (0.940, 0.970) GeV/c 2 , respectively.

A. Single-tag yields
The ST D signals are identified using the beam-constrained mass, where − → p Dtag is the momentum of the D candidate. To improve the signal purity, the energy difference ∆E = √ s/2 − E Dtag for each candidate is required to be within approximately ±3σ ∆E around the ∆E peak, where σ ∆E is the ∆E resolution and E Dtag is the reconstructed ST D energy. The explicit ∆E requirements for all reconstructed ST modes are listed in the second column of Table II. For the ST channels of K + π − , K + K − and π + π − , backgrounds of cosmic rays and Bhabha events are removed with the following requirements. First, the two charged tracks must have a TOF time difference of less than 5 ns and they must not be consistent with being a muon pair or an e + e − pair. Second, there must be at least one EMC shower with an energy larger than 50 MeV or at least one additional charged track detected in the MDC.
The M BC distributions for the ST modes are shown in Fig. 2. To obtain the ST yields reconstructed by these modes, maximum likelihood fits are performed to these spectra, where the signal peak is described by a Monte Carlo (MC) simulated shape convolved with a double-Gaussian function, and the combinatorial background is modeled with an ARGUS function [35]. In addition to the combinatorial background, there are also some peaking backgrounds in the signal region of M BC . These peaking backgrounds are included in the yields obtained from fits to M BC spectra and hence must be subtracted. For example, for the ST modes of K + π − , K + π − π 0 and K + π − π − π + , there are small contributions of wrong-sign (WS) peaking backgrounds in the STD 0 samples, which originate from the DCS-dominated decays of D 0 → K + π − , K + π − π 0 and K + π − π − π + . In addition, the D 0 → K 0 S K + π − (K 0 S → π + π − ) decay is a source of WS peaking background for the ST decayD 0 → K + π − π − π + . Overall, the peaking background contamination rates are less than 1% for the ST modes of K + π − , K + π − π 0 and K + π − π − π + . For the CP -eigenstate ST channels K 0 S π 0 (π 0 ) and π + π − π 0 , the peaking background rates are 0.8%(3.9%) and 3.9%, dominated by the D meson decays to π + π − π 0 (π 0 ) and K 0 S π 0 , respectively. The D → K 0 S π + π − π 0 decay forms the dominant peaking backgrounds and accounts for contamination rates of 13.7%, 6.3% and 3.8% in the fitted ST yields for K 0 S ω, K 0 S η π + π − π 0 and K 0 S η γπ + π − , respectively. Additionally, the sample of ST K 0 S π + π − decays includes a 2% contamination from the peaking-background D → π + π − π + π − . The sizes of these peaking backgrounds are all estimated from MC simulation and then subtracted from the fitted ST yields. The background-subtracted yield and the efficiency for each of the ST modes are summarized in the third and fourth columns of Table II, respectively. The ST efficiencies are determined from the simulated data where one D meson is forced to decay to the reconstructed final states and the other D meson is allowed to decay to any final state. The values of ST vary from ∼65% for decay modes with two charged particles in the final state to ∼13% for final states with multiple composite and neutral particles such as K 0 S η π + π − η . The ST yields of the modes K + e −ν e , K 0 L π 0 and K 0 L π 0 π 0 , which cannot be directly reconstructed, are estimated from knowledge of the number of neutral D meson pairs N DD , the estimated ST efficiencies ST tag , and their BFs B tag reported in Ref. [33], where the D → K 0 S π 0 π 0 BF is used as a proxy for D → K 0 L π 0 π 0 . The yields are calculated from the relations . The ST efficiencies, ST tag , of detecting these three decays are estimated by evaluating the ratios between the corresponding DT (discussed later in Sec. IV F) and ST efficiencies, which are determined to be 61.35%, 48.97% and 26.20% for D → K + e −ν e , D → K 0 L π 0 and D → K 0 L π 0 π 0 , respectively. The ST yields of D → K − e + ν e , D → K 0 L π 0 and D → K 0 L π 0 π 0 are also included in Table II, in which the uncertainties from the BFs, N DD and the detection efficiencies are presented.
In those cases where the decay products of the tag mode are fully reconstructed and the signal mode is D → K 0 S π + π − , the signal decay is built by using the other tracks in the event recoiling against the ST D meson. The same selection on track parameters and the K 0 S candidate is imposed as de- candidate, is required to be between −30 and 33 MeV. If multiple combinations are selected, the one with the minimum |∆E | is retained. The beam-constrained mass is defined as M sig where p sig is the momentum of the signal-decay candidate.
The DT yield is determined by performing a twodimensional unbinned maximum-likelihood fit to the M sig BC (signal) vs. M tag BC (tag) distribution. An example distribution for the tag mode D → K + π − is shown in Fig. 3. The signal shape of the M sig BC vs. M tag BC distributions is modeled with a two-dimensional shape derived from simulated data convolved with two independent Gaussian functions representing the resolution differences between data and simulation. The parameters of the Gaussian functions are fixed at the values obtained from the one-dimensional fits of the M sig BC and M tag BC distributions in data, respectively. The combinatorial backgrounds in the M sig BC and M tag BC distributions are modeled by an ARGUS function in each dimension where the parameters are determined in the fit. The events that are observed along the diagonal arise from mis-reconstructed DD decays and from qq events. They are described with a product of a double-Gaussian function and an ARGUS function rotated by 45 • [36]. The kinematic limit and exponent parameters of the rotated ARGUS function are fixed, while the slope parameter is determined by the fit. The peaking backgrounds in the M sig BC and M tag BC distributions are described by using a shape derived from simulation convolved with the same Gaussian function as used for the signal. The decay D → π + π − π + π − , which accounts for about 2% peaking background to D → K 0 S π + π − signal, is predominantly CPeven [37], and hence the yields of this peaking background are adjusted from the expectation of simulation to account for the effects of quantum correlation. Figure 4 shows the projections of the two-dimensional fits on the M sig BC distribution for all the fully reconstructed ST decay modes.
The DT yield of K 0 S π + π − vs. K 0 S π + π − is crucial for determining the s i values and thus it is desirable to increase the reconstruction efficiency for these events. Therefore three independent selections are introduced in order to maximize the yield of D → K 0 S π + π − vs. D → K 0 S π + π − candidates. The first selection requires that both K 0 S π + π − final states on the signal and tag side are fully reconstructed. However, in order to increase the efficiency, the PID requirements on the pions originating from both the signal and tag D mesons are removed and the K 0 S candidate needs only satisfy L/σ L > 0. This looser selection is applied to both D mesons and allows for an increase in yield of approximately 20% with only a slight increase in background.
The second selection class allows for one pion originating from the D meson to be unreconstructed in the MDC, denoted as K 0 S π + π − miss . Events with only three remaining charged tracks recoiling against the D → K 0 S π + π − ST are searched for. The K 0 S and pion are identified with the same criteria used to select the ST candidates. The missing pion is inferred by calculating the missing-mass squared (M 2 miss ) of the event, which is defined as where p sig is the momentum of the fully reconstructed D → K 0 S π + π − candidate and i E i and i p i are the sum of the energy and momentum of the other reconstructed particles that form the partially reconstructed D meson candidate. Throughout this paper, in order to determine the signal yields of the DT containing a missing particle, an unbinned maximum-likelihood fit is performed to the defined kinematic distribution, i.e. M 2 miss (or U miss discussed in Sec. IV D). The signal and background components are described using shapes from simulated data where the signal shape is further convolved with a Gaussian function. The relative yields of the peaking backgrounds to the signals are fixed in the fits from information of the simulated data. Figure 5(a) shows the M 2 miss distribution from the partially reconstructed D → K 0 S π + π − vs. D → K 0 S π + π − miss candidates. The distribution peaks at M 2 miss ∼ 0.02 GeV 2 /c 4 , which is consistent with the missing particle being a π ± . The peaking backgrounds are approximately 3% of the signal yield and are primarily from the D → π + π − π + π − decay.
The third D → K 0 S π + π − vs. D → K 0 S π + π − selection identifies those events where one K 0 S meson decays to a π 0 π 0 pair. Events where there are only two remaining oppositelycharged tracks recoiling against the ST D → K 0 S π + π − are selected and these tracks are classified as the π + and π − from the D meson. To avoid the reduced efficiency associated with reconstructing both π 0 mesons from the K 0 S , only one of the them is searched for. This type of tag is referred to as K 0 S (π 0 π 0 miss )π + π − . The missing-mass squared of the event is defined in the same way as in Eq. (13) and the summation is over the π + , π − , and π 0 mesons that are reconstructed on the tag side. A further variable, M 2 miss , where the reconstructed π 0 is also not included in the summed energies and momenta of the tag-side particles is also computed. For true D → K 0 S π + π − decays this variable should be consistent with the square of the K 0 S meson nominal mass. Therefore, candidates that do not satisfy 0.22 < M 2 miss < 0.27 GeV 2 /c 4 are removed from the analysis in order to suppress background from D → π + π − π 0 π 0 decays. Figure 5  sultant M 2 miss distribution of the accepted candidates in data. There remains a contribution of peaking background dominated from D → π + π − π 0 π 0 decays, where the rate relative to signal is determined from simulated data to be around 15%.
DT candidates are also reconstructed with the missing-mass squared technique as the K 0 L particle is not directly detectable in . The listed uncertainties are statistical only.

Mode
ST DT  the BESIII detector. In the rest of the event containing a D → K 0 S π + π − ST, a further π 0 or π 0 π 0 pair is reconstruct- candidates, respectively. A peak at the square of the mass of the K 0 L meson is clearly visible. In this case the peaking backgrounds come from events where the decay products of the K 0 S have not been reconstructed and therefore the K 0 S meson has been identified as a K 0 L meson. The peaking backgrounds from D → K 0 S π 0 and D → K 0 S π 0 π 0 comprise 5% and 9%, respectively, of the signal sample.
The D 0 → K − e + ν e vs.D 0 → K 0 S π + π − DT candidates are reconstructed by combining an ST K 0 S π + π − candidate with a K − and a positron candidate from the remaining tracks in the event. Events with more than two additional charged tracks that have not been used in the ST selection are vetoed. Information concerning the undetected neutrino is obtained through the kinematic variable where E K and E e are the energy of the kaon and electron from the semi-leptonic D decay candidate, and p miss is the missing momentum carried by the neutrino. The momentum 500 1000 -π + π 0 S K vs. p miss is defined as p miss = p sig − p K − p e . Figure 5(e) shows the U miss distribution for D 0 → K − e + ν e candidates in data, where a peak centered on U miss = 0 is observed due to the negligible mass of the neutrino.
E. Double tags with K 0 L π + π − To identify the signal candidates from D → K 0 L π + π − decays, only two additional and oppositely charged good tracks are required in an event where one of the ST has been selected. These two tracks are identified as the π + and π − from the D meson. Events that contain any additional charged tracks with the distance of closest approach to the IP less than 20 cm along the beam direction are vetoed. This requirement reduces background from K 0 S → π + π − decays. To reject the backgrounds containing π 0 and η mesons, events are vetoed where the invariant mass of any further photon pairs are within the ranges (0.098, 0.165) GeV/c 2 and (0.48, 0.58) GeV/c 2 . This requirement retains about 80% of the signal while reducing more than 90% of the peaking backgrounds from D → K 0 S π + π − , where K 0 S → π 0 π 0 . The residual peaking background rate in D → K 0 L π + π − selected candidates is 5% of the signal yield and is primarily from the decay D → K 0 S (π 0 π 0 )π + π − . Figure 6 shows the M 2 miss distributions of the accepted D → K 0 L π + π − candidates in data.

F. Dalitz plot distributions
The DT yields of K 0 S π + π − and K 0 L π + π − tagged by different channels are shown in the fifth and seventh columns of Table II, respectively. Their selection efficiencies ( DT ) are also listed in the sixth and eighth columns of Table II. The DT selection efficiencies are determined in simulation where the signal and tag D meson are both forced to decay to the final states in which they are reconstructed. The efficiency is determined as the number of DT candidates selected divided by the number of events generated.
The DT yields of D → K 0 S(L) π + π − involving a CP eigenstate are a factor of 5.3(9.2) larger than those report-ed in Ref. [23]. The yields of K 0 S π + π − tagged with D → K 0 S(L) π + π − decays are a factor of 3.9(3.0) larger than those in Ref. [23]. These increases come not only from the larger data set available at BESIII but also from the additional tag decay modes and partial reconstruction selection techniques.
The resolutions of M 2 K 0 S π ± and M 2 K 0 L π ± on the Dalitz plot are improved by requiring that the two neutral D mesons conserve energy and momentum in the center-of-mass frame, and the decay products from each D meson are constrained to the nominal D 0 mass [33]. In addition the K 0 S decay products are constrained to the K 0 S nominal mass [33]. Finally, the missing mass of K 0 L candidates is constrained to the nominal value [33]. The study of simulated data indicates that the resulting resolutions of M 2 respectively. It should be noted that the finite detector resolution can cause the selected events to migrate between Dalitz plot bins after reconstruction, which should be incorporated in evaluating the expected DT candidates observed in Dalitz plot bins. More details are presented in Secs. V B and V C.
The Dalitz plots for D 0 → K 0 S π + π − and D 0 → K 0 L π + π − vs. the flavor tags selected from the data are shown in Fig. 7. In order to merge the D 0 andD 0 decays the exchange of coordinates M 2 Figure 7 also shows the CP -even and CP -odd tagged signal channels selected in the data. The effect of the quantum correlation in the data is immediately obvious by studying the differences in these plots. Most noticeably, the CP -odd com-  The fit used to determine the strong-phase parameters is based on the Poisson probability to observe N events in a S π + π − and K 0 L π + π − events in data.
phase space region given the expectation value N . To measure the observed yields, the data are divided into the phase space regions based on their Dalitz plot coordinates (m 2 + , m 2 − ). A small fraction of candidates (∼0.3%) fall outside the defined bins. This is because the knowledge of the D 0 mass has improved since the model used to define the phase space regions was determined. This improvement leads to a slightly larger allowed phase space in the current analysis compared to the maps of the phase space regions. These outlying candidates are assigned to the bins to which they are closest.
In the K 0 S,L π + π − Dalitz plots of the flavor-tagged samples, the positive and negative bins are distinguishable, and hence yields are measured in 16 bins for each final state. In contrast, the CP -tagged Dalitz plots are symmetric about the line m 2 + = m 2 − (see Eqs. (7) and (10)) and so the entries are summed for bins i and −i. Exploiting this symmetry reduces the statistical fluctuations for those CP tags where the yields are low.
The K 0 S(L) π + π − vs. K 0 S π + π − samples are described by two Dalitz plots. Therefore it is necessary to determine the yields for the ith bin (D i ) of one plot and the jth bin (D j ) of the other, in order to obtain the quantities M ij (M ij ) that occur in Eq. (9)  It follows that events can be classified into those where both decays occur in the Dalitz plots on the same side of the m 2 + = m 2 − line, and those when they are on different sides. For the case where both D mesons are fully reconstructed as D → K 0 S π + π − , it is not possible to distinguish between D i and D j , and thus M ij is combined with M ji . The partially reconstructed D → K 0 S π + π − samples are treated in the same way, despite the distinguishability of the final states, in order to avoid low yields. In the K 0 L π + π − vs. K 0 S π + π − sample D i is chosen to specify the K 0 S π + π − Dalitz plot bin, and D j the K 0 L π + π − bin. In this case M ij and M ji are distinguishable and cannot be combined. Following these considerations, the samples with two Dalitz plots are divided into 72 and 128 bins.
In each bin of phase space there are candidates that are from signal, combinatorial background, and peaking backgrounds. The yields for each DT mode are determined in the same way as in Sec. IV, although in some regions where the yields are low it is necessary to fix some parameters from the fit to data over the full phase space. The observed combinatorial background yield is determined in the fit and not considered further. Although the expected peaking-background yield can be calculated with MC simulation, the fit cannot distinguish the observed peaking background yield from the signal yield. Therefore the observed yield N obs in each phase space region is the sum of signal and peaking background.

B. Determination of Ki and K i
The yields of K i and K i are necessary to determine the expected yields in the decays sensitive to the strong-phase parameters. As discussed in Sec. IV F, the finite detector resolution can cause the individual decays to migrate between Dalitz plot bins after reconstruction. Furthermore, the migration effects between D 0 → K 0 S π + π − and D 0 → K 0 L π + π − are also different due to their resolution differences. Studies indicate that neglecting bin migration induces average biases of 0.7 (0.3) times the statistical uncertainty in the determination of c i (s i ), and hence it is important to correct for this effect in the analysis.
To account for this effect, the number of observed signal events (N obs( ) i ) for flavor-tagged K 0 S(L) π + π − decays in the ith bin of the Dalitz plot is written where ij is the efficiency matrix which describes the reconstruction efficiency and migration effects across Dalitz plot bins associated with reconstruction of tag and signal decays. The efficiency matrix ij can be obtained by analyzing a sample of signal MC events which are generated as e + e − → ψ(3770) → D 0D0 , where theD 0 meson decays to the ST modes and D 0 → K 0 S(L) π + π − . The efficiency ma-trix ij for detecting D → K 0 S(L) π + π − decay is given by where N rec ij is the number of signal MC events generated in the jth Dalitz plot bin and reconstructed in the ith Dalitz plot bin, N gen j is the number of signal MC events which are generated in the jth Dalitz plot bin and ST is the ST efficiency. An example of the efficiency matrix ij for K 0 S(L) π + π − vs. K + π − in the equal ∆δ D binning scheme is shown in Table III. Thus, the value of K ( ) i in the ith Dalitz plot bin for D 0 → K 0 S(L) π + π − decay is obtained by In addition, the migration effects in the ith Dalitz plot bin can be estimated by using R i = ii / j ij , which denotes the fraction of the reconstructed events falling outside the true Dalitz plot bins. From the efficiency matrix ij listed in Table III, it is estimated that the bin migration effects range within (3-12)% and (3-18)% for the K 0 S π + π − and K 0 L π + π − signals with the equal ∆δ D binning scheme, respectively.
Moreover, the event yields of K 0 S π + π − and K 0 L π + π − selected against hadronic flavored tags are also contaminated by DCS decays [22,23]. To account for this effect, the flavortagged yield in each Dalitz plot bin is scaled by a correction factor f ( ) i (f i for K 0 S π + π − and f i for K 0 L π + π − ). The correction factors for the hadronic tags K + π − , K + π − π 0 and K + π − π − π + are calculated by Here r F D is the ratio of the DCS amplitude to the Cabibbofavored decay amplitude and δ F D is the corresponding strongphase difference. For multibody final states these two quantities are averaged over the decay phase space. The coherence factor [38] R F equals unity for two-body decays, and has been measured for D → K + π − π − π + and D → K + π − π 0 [39,40]. The values of these parameters and the corresponding references are given in Table IV. Furthermore, f (m 2 + , m 2 − ) and f (m 2 + , m 2 − ) are the amplitudes of D 0 → K 0 S π + π − and D 0 → K 0 L π + π − , respectively. The amplitude of the decay D 0 → K 0 S π + π − is taken from the model given in Ref. [42]. The decay amplitude of D 0 → K 0 L π + π − has not been studied, however a decay model can be estimated by adjusting the D 0 → K 0 S π + π − model in the same way as dis-cussed in Refs. [22,23], and using that K 0 S and K 0 L mesons are of opposite CP , to an excellent approximation. Starting with the D 0 → K 0 S π + π − model in Ref. [42], the Cabibbofavored amplitudes are unchanged and the amplitudes of the DCS components gain a factor −1. For the CP -eigenstate amplitudes, such as K 0 S ρ(770) 0 , the D 0 → K 0 L π + π − amplitude can be related to the D 0 → K 0 S π + π − amplitude by multiplying the latter by a factor (1 − 2re iδ ) [43], where r is of the order of tan 2 θ C . Here, θ C is the Cabibbo angle and δ is an unknown phase. To determine central values for f i , the parameters r and δ are varied a number of times where δ is assumed to have an equal probability to lie between 0 • and 360 • and r is assumed to have a Gaussian distribution with mean tan 2 θ C and width 0.5 × tan 2 θ C . The mean value of TABLE III. Efficiency matrix ij (%) for K 0 S,L π + π − vs. K + π − in the equal ∆δD binning scheme. The column gives the true bins j, while the row gives the reconstructed bin i, in which the decay BF of K 0 S → π + π − is not included.
is computed for each flavor-tag. Figure 8 shows the measured values of F i and F i in various Dalitz plot bins for the flavor-tagged D 0 → K 0 S π + π − and D 0 → K 0 L π + π − events, respectively, observed in data. The measured values of F i and F i are consistent between the different categories of flavor tags, which provides a good validation for the extracted K i and K i in data. In order to recover the summed K i from all flavor-tagged K 0 S π + π − with K 0 S → π + π − used in this analysis, the values of F i in Table V should be multiplied by 58607, 58647, and 58595 for the equal ∆δ D , optimal, and modified binning schemes, respectively. To obtain the values of K i the corresponding factors are 80718, 80661, and 80706.

C. Expected DT yields in Dalitz plot bins
The expected yields, N , are a sum of expected signal and peaking background contributions. The expected signal yields are calculated from the equations given in Sec. II, with adjustments made to account for bin migration and selection and reconstruction efficiencies, so that they can be compared to (left) F i/−i for K 0 S π + π − and (right) F i/−i for K 0 L π + π − from (top) equal ∆δD, (middle) optimal and (bottom) modified optimal binning schemes, respectively. The horizontal and vertical error bars denote the bin intervals and the statistical uncertainties of the event fractions, respectively. the yields from data. For CP -tagged D → K 0 S π + π − decays, the expected signal yield in the ith bin is given by where K 0 S π + π − ij is the efficiency matrix for detecting D → K 0 S π + π − vs. the particular CP tag under consideration defined similarly as in Eq. (17). The efficiency matrix is defined to take into account the merging of the ith and −ith regions in data (i.e. it has size 8 × 8). The normalization factor h CP is defined as S CP /S FT . For CP -tagged D → K 0 L π + π − , the expected signal yield in the ith bin is given by where h CP is given by S CP /S FT , and S FT is the sum of ST yields for the three hadronic flavor tags used to determine the values of K i .
For the DT D → K 0 S π + π − vs. D → K 0 S π + π − the efficiency matrix, , is a 72 × 72 matrix where each value of the indices n and m corresponds to one of the 72 distinct bin-pairs. The expected signal yields are expressed where i m and j m correspond to the ith and j th bins of the mth bin-pair and h corr = (N DD /S 2 FT ) × αβ. The constant α results from the symmetry relations used to combine bin pairs. The value of α is 1 when the i, j values of the index n satisfy |i| = |j| and 2 otherwise. The constant β arises from the symmetry of the signal and tag decays and has value 1 for the selections where both K 0 S mesons decay via K 0 S → π + π − , and value 2 × B(K 0 S → π 0 π 0 )/B(K 0 S → π + π − ) when one K 0 S meson decays to the π 0 π 0 final state. Here, B(K 0 S → π 0 π 0 ) and B(K 0 S → π + π − ) are BFs for K 0 S → π 0 π 0 and K 0 S → π + π − , respectively. For the DT D → K 0 L π + π − vs. D → K 0 S π + π − the expected signal yields are expressed as where h corr = 2N DD /(S FT S FT ). The DT efficiency matrix of detecting D → K 0 L π + π − vs. D → K 0 S π + π − , K 0 L π + π − vs.K 0 S π + π − nm , is a 128 × 128 matrix where each value of the indices n and m corresponds to one of the 128 distinct bin-pairs.
The peaking-background yields integrated over the Dalitz plot have been estimated for each DT in Sec. IV. The majority of the peaking backgrounds are CP eigenstates, for example the decay D → K 0 S π 0 forms a peaking background to the tag D → K 0 L π 0 . The simulated data cannot accurately describe the distribution of the peaking background over the Dalitz plot since it does not account for quantum correlations. For the peaking backgrounds which are CP eigenstates the expected yields are distributed according to Eq. (7) with an appropriate normalization factor to take into account the expected yield integrated over the Dalitz plot. The values of c i and s i used to make the initial estimate of the expected peaking background yields are computed from the D → K 0 S π + π − amplitude model [42]. As the yields of peaking backgrounds are small compared to the signal yields the effects of migration or small variations in efficiency over the Dalitz plot are ignored. A similar estimate for the peaking background of D → K 0 S π + π − vs. D → K 0 S π + π − in the D → K 0 L π + π − vs. D → K 0 S π + π − DT can be estimated through Eq. (9). The peaking background D → K 0 S π + π − π 0 in the D → K 0 S ω tag is treated as CP -odd, as indicated by the results in Ref. [44]. The strong-phase parameters of D → π + π − π 0 π 0 are not known and in this case the peaking background is distributed as observed in simulated data. The expected distribution of the remaining peaking backgrounds that occur at low rates, such as D → π + π − π + π − , are also taken from simulation.
where P (N obs , N exp ) is the Poisson probability to observe N obs events given the expected number N exp . The observed yield of signal and peaking background in the ith bin or nth bin-pair is denoted N obs i,n , and N exp i,n is defined to account for both expected signal and peaking background from the same region. Biases can occur in the case where N obs i,n is close to zero. To mitigate this effect the observed and expected yields of the three selections of D → K 0 S π + π − vs. D → K 0 S π + π − DT candidates are summed together. The observed and expected yields of the two final states of the D → K 0 S η tag are also added together and the same is done for both final states of the D → K 0 S η tag. The χ 2 term in Eq. (24) is which constrains the measured differences c i − c i (s i − s i ) to the predicted differences, ∆c i (∆s i ), where δ∆c i (δ∆s i ) are the uncertainties in the predictions. The presence of the constraint is necessary in order to improve the precision of s i and s i , and introduces very weak model assumptions in the fit. The expected values of c i and s i are determined from the D → K 0 S π + π − amplitude model in Ref. [42]. The expected values of c i and s i are determined from the assumed D → K 0 L π + π − amplitude model described in Sec. V B, where the central values come from the mean of the strongphase distributions generated using different values of r and δ. In order to determine δ∆c i and δ∆s i , the values of ∆c i and ∆s i are also estimated using the models of D → K 0 S π + π − reported in Refs. [45,46] with the same transformation to estimate the D → K 0 L π + π − decay model. In order to assign δ∆c i and δ∆s i , the larger deviation of the central values of ∆c i and ∆s i using these two alternative models is taken as part of the uncertainty and added in quadrature to the uncertainty from the choice of r and δ. The CP -tagged data can also be used to fit only c i and c i where the likelihood does not contain a constraint on the difference between these parameters. The measured differences from this fit are consistent with the predicted values of ∆c i , which gives further confidence in the transformations used to define the D → K 0 L π + π − amplitude model. Table VI summarizes the expected c i , s i , ∆c i and ∆s i and uncertainties for the three binning schemes.
To resolve the ambiguity in the sign of s i present in Eqs. (22) and (23), the starting values of the parameters of the fit are set to be consistent with the model prediction. An iterative fit is performed to the data. After each iteration the expectation values of the peaking backgrounds that use c i and s i as input are recalculated using the c i and s i values determined by the fit. Three iterations are required to provide a stable result. The fitted strong-phase parameters c i , s i , c i , and s i are summarized in Table VII Furthermore, several checks are performed to assess the stability of the fit results. The fits are repeated on different subsets of the data, for example, separating partially and fullyreconstructed K 0 S π + π − events. Further tests involve removing specific tags, such as π + π − π 0 , K 0 L π 0 and K 0 L π 0 π 0 . The results from these checks are consistent with the default values of c ( ) i and s ( ) i . Furthermore, the c i and s i results are found to be robust when K 0 L π + π − tags are removed from the fit.

VI. SYSTEMATIC UNCERTAINTIES
Uncertainties associated with the selection, tracking and PID efficiencies do not bias the measurement as the expected DT yields are calculated using the ST yields and determined values of K ( ) i . This use of data-driven quantities to provide the normalization means that detector effects on the common selection affect the observed and expected DT yields in the same way, and hence these systematic uncertainties are not considered further.
The statistical uncertainties on the measured K  Table IV and by assessing the impact of using the alternative D → K 0 S π + π − models reported in Refs. [45] and [46] to calculate the f A difference in the resolution between simulation and data introduces an uncertainty in the efficiency matrices. The dif-ference in resolution is quantified by studying the mass spectrum of the K * (892) resonance found in D → K 0 S π + π − decays. The mass spectrum is fitted with a shape determined by simulation convolved with a Gaussian function, which defines the difference in resolution between data and simulation. The Gaussian has a mean of 0.23 MeV/c 2 and width 0.21 MeV/c 2 . The variables M K 0 S π − and M K 0 S π + of all simulated events used to determine the efficiency matrices are smeared by a Gaussian with these parameters and new efficiency matrices are calculated. The same procedure is performed on the mass spectrum of K * (892) from D → K 0 L π + π − decays and the differences here are described with a Gaussian with mean 4.0 MeV/c 2 and width 2.0 MeV/c 2 . The fit to determine the strong-phase parameters is repeated with the new efficiency matrices and the differences between these fit results and the nominal values are assigned as the systematic uncertainty due to residual differences between the momentum resolution in data and simulation. The impact of finite samples of simulated data to determine the efficiency matrices on the strong phases is assessed by varying the matrix elements by their statistical uncertainties. This is repeated multiple times, and the data are refitted using these new matrices to determine the expected yields. The resulting width of the distribution of the values of the strong-phase parameters is assigned as the systematic uncertainty due to the size of the simulated samples.
The expectation values of the peaking background have systematic uncertainties due to the inputs used to calculate their integrated yields and the assumptions concerning the distribution over the Dalitz plot. For the uncertainty from the integrated yields, the expected yield of peaking background in each phase space region is varied according to a Gaussian distribution. This distribution has the nominal value of the peaking-background yield as the mean, and a width which combines the uncertainties from the BFs of peakingbackground decays, and the uncertainties arising from tracking [47], PID [47], and π 0 reconstruction efficiencies [48]. The distributions over the Dalitz plots for peaking backgrounds that are CP eigenstates, or D → K 0 S π + π − for D → K 0 L π + π − signals, are dependent on the values of c i and s i . As the iterative fit procedure recalculates the peaking background with updated values of c i and s i , no further systematic uncertainty is assigned for these backgrounds. The D → π + π − π 0 π 0 peaking background constitutes a significant contribution to the observed yields in D → K 0 S π + π − vs. D → K 0 S π + π − where one K 0 S meson decays to the π 0 π 0 final state. To find an alternative distribution of this background a DT sample of D → K 0 S π + π − vs. D → π + π − π 0 π 0 events is fully reconstructed in data. The distribution over the Dalitz plot is found by assigning the K 0 S mass to the π 0 pair. This distribution is used instead of the nominal one (from simulation) in the fit and small shifts are observed in the strongphase parameters that are assigned as an additional contribution to the systematic uncertainties arising from the DT peaking backgrounds.
The effects from D 0D0 mixing are not considered in the nominal fit. The required correction factor for CP ± eigenstate ST yields is 1/(1 − η ± y D ), where η ± = ±1 and the mixing parameter y D = (0.62±0.08)% [33]. The data are fitted using the corrected ST yields and the difference with respect to the nominal results is assigned as the systematic uncertainty due to charm mixing.
Systematic uncertainties in the observed DT yields arise from the fit procedure and the description of combinatorial background. Due to the low candidate yields in multiple of the phase space regions, small biases in the fitted yields can be present. The sizes of these biases are determined in pseudoexperiments. An alternative combinatorial background shape is employed and the difference in N obs between this fit and the nominal is added in quadrature to the bias estimate to determine the systematic uncertainty on the observed yields. All the observed yields are smeared within these uncertainties and the fit is repeated. The resulting width of the distribution of values of the strong-phase parameters is assigned as the systematic uncertainty due to the DT yields.
The systematic uncertainties of the measured strong-phase parameters c i , s i , c i and s i for the equal ∆δ D , optimal, and modified optimal binning schemes are summarized in Tables VIII, IX and X, respectively. There is no source of systematic uncertainty that is dominant for all strong-phase parameters. The statistical uncertainty obtained from the fit includes the contribution related to the associated uncertainties on ∆c i and ∆s i through the χ 2 term of Eq. (24). In order to estimate this contribution, the fit is repeated in a configuration where ∆c i and ∆s i are fixed. The difference in quadrature between the uncertainties from this fit and the nominal approach provides an estimate of the contribution to the uncertainty from the constraint. This estimate is also given in Tables VIII, IX and X, and it is seen that this contribution to the overall uncertainty is small. The measurements of the strong-phase parameters are limited by their statistical uncertainties. The correlation matrices for the statistical and systematic uncertainties associated with different binning schemes are given in Tables XI−XVI. The measurements are displayed in Fig. 9, together with the model predictions from Ref. [42], which are seen to be in reasonable agreement. Given the compatibility between the current measurements and those reported by the CLEO collaboration [23] an additional set of fits is performed, where the CLEO results are imposed as a Gaussian constraint in Eq. (24). These results are presented in Appendix A.

VII. IMPACT ON γ/φ3 MEASUREMENT
The model-independent measurement of γ described in Ref. [4] is performed by comparing the number of B − → DK − , D → K 0 S π + π − events in a given Dalitz plot bin with the integral of the square of the amplitude given in Eq. (1) over the same region. An analogous expression for the B + events is also used. Therefore the expected yield of B ± events in a Dalitz plot region is a function of K i , c i , s i , and γ, δ B and r B , the underlying parameters of interest, and is given by i measured in this work (red dots with error bars), the expected values from Ref. [42] (blue open circles) as well as CLEO results (green open squares with error bars) in Ref. [23]. The top plots are from the equal ∆δD binning, the middle plots from the optimal binning and plots from the modified optimal binning scheme are on the bottom. The circle indicates the boundary of the physical region c In order to assess the impact of the uncertainty in the strongphase parameters on a measurement of γ, a large simulated data set of B ± events is generated according to the expected distribution given the measured central values of K i , c i , and s i and the input values γ = 73.5 • , r B = 0.103, and δ B = 136.9 • , which are close to the current central values of these parameters from existing measurements [41]. The simulated data are fit many times to determine γ, δ B and r B . The values of c i and s i used in each fit are sampled from the measured values smeared by their uncertainties, where the correlations between the measurements are taken into account. The uncertainty on the measured K i is not considered since experiments are expected to use their own data to provide this input [13]. The overall yield of the generated B ± sample is sufficiently large to ensure that the statistical uncertainty from the fit is negligible. Therefore the width of the distribution of the fitted value of γ is an estimate of the uncertainty on γ due to the precision of the strong-phase parameters. The distribution of the fitted value of γ in the three binning schemes is shown in Fig. 10.
Based on this study, the uncertainty on γ due to the measured uncertainty on c i and s i is found to be 0.7 • , 1.2 • and 0.8 • for the equal ∆δ D , optimal and modified optimal binning schemes, respectively. The very small phase space regions in the optimal binning scheme are the cause for the larger propagated uncertainty in this case. Very small biases of less than 0.2 • are observed due to some values in the fit being unphysical, i.e. c 2 i + s 2 i > 1. The size of the uncertainty on γ is approximately a factor of three smaller than from the CLEO measurements [23]. The predicted statistical uncertainties on γ from LHCb prior to the start of High-Luminosity LHC operation in the mid 2020s, and from Belle II is expected to be 1.5 • [49,50]. Therefore the uncertainty associated to the strong-phase measurements presented here will not be dominant in the determination of γ for Belle II or for LHCb until then. The measurements of c i and s i can also be used for determination of strong-phase parameters in other multi-body decay modes of D mesons, where the D 0 → K 0 S π + π − decay is used as a tag [26,37,39,44,51]. Here, the improved precision leads to smaller systematic uncertainties on the strongphase parameters in other D decay modes, which subsequently reduces associated systematic uncertainties on γ when these D decay modes are use to measure γ in B ± → DK ± decays.

VIII. SUMMARY
Measurements of the relative strong-phase differences between D 0 andD 0 → K 0 S,L π + π − in bins of phase space have been performed using 2.93 fb −1 of data collected at √ s=3.773 GeV collected with the BESIII detector. These results are on average a factor of 2.5 (1.9) more precise for c i (s i ) and a factor of 2.8 (2.2) more precise for c i (s i ) than the previous measurements of these parameters [23]. This improvement arises from the combination of a larger data sample, an increased variety of CP tags, and broader use of the partial reconstruction technique to improve efficiency. The strongphase parameters provide an important input in a wide range of CP violation measurements in the beauty and charm sectors. The propagated uncertainty from these measurements on the CKM parameter γ determined through the analysis of B ± → D K 0 S π + π − K ± events is expected to be 0.7 • , 1.2 • and 0.8 • for the equal ∆δ D , optimal and modified optimal binning schemes, respectively. This improved precision will ensure that measurements of γ from LHCb and Belle II over the next decade are not limited by the knowledge of these strong-phase parameters, and also be invaluable in studies of charm mixing and CP violation.  As the results presented here and those from the CLEO collaboration [23] are compatible it is legitimate to combine them in order to provide a single set of results that benefits from both measurements. The combination is performed by performing the fit described in Sec. V to the double tags with an additional multi-dimensional Gaussian constraint present on the strong-phase parameters. This constraint comes from the central values and the covariance matrices in Ref. [23].
A small, additional contribution to these covariance matrices, determined through pseudo-experiments, accounts for the effects reported in Ref. [40].
The systematic uncertainties reported in Tables VIII − X are added in quadrature to those from the fit, which include contributions from the BESIII statistical and CLEO statistical and systematic uncertainties. The central values and their uncertainties for the three binning schemes are reported in Table XVII and Tables XVIII, XIX, and XX show the associated covariance matrices.