Improved model-independent determination of the strong-phase difference between $D^{0}$ and $\bar{D}^{0}\to K^{0}_{\mathrm{S,L}}K^{+}K^{-}$ decays

We present a measurement of the strong-phase difference between $D^{0}$ and $\bar{D}^{0}\to K^{0}_{\rm S,L}K^{+}K^{-}$ decays, performed through a study of quantum-entangled pairs of charm mesons. The measurement exploits a data sample equivalent to an integrated luminosity of 2.93~fb$^{-1}$, collected by the BESIII detector in $e^{+}e^{-}$ collisions corresponding to the mass of the $\psi(3770)$ resonance. The strong-phase difference is an essential input to the determination of the Cabibbo-Kobayashi-Maskawa (CKM) angle $\gamma/\phi_3$ through the decay $B^{-}\to DK^{-}$, where $D$ can be either a $D^{0}$ or a $\bar{D}^{0}$ decaying to $K^{0}_{\rm S,L}K^{+}K^{-}$. This is the most precise measurement to date of the strong-phase difference in these decays.


I. INTRODUCTION
In the Standard Model (SM), the charged-weak interaction in the quark sector is described by the Cabibbo-Kobayashi-Maskawa matrix (CKM) [1]. One of the primary goals of flavor physics experiments is to determine the angles α, β and γ (or φ 2 , φ 1 and φ 3 ) of the b − d CKM unitary triangle precisely. Currently, the most precise measurements of γ are extracted using tree-level B − → DK − decays [2]. Here and elsewhere in this paper D refers to either a D 0 or aD 0 meson decaying into the same final state and charge conjugation is implicit, unless stated otherwise. The sensitivity to γ arises from the interference of two amplitudes: b → cūs that results in the B − → D 0 K − decay, and b → ucs that leads to the B − →D 0 K − decay. The latter amplitude is both CKM-and color-suppressed relative to the former. The value of γ measured with such tree-level transitions is insensitive to loop-level contributions [3]. Therefore, tests for new physics that are made by comparing unitarity triangle parameters measured using tree and loop processes can be improved by more precise determinations of γ [4,5].
Different methods of determining γ are classified based upon the decay products of the D decay: CP eigenstates (GLW method) [2], flavor-eigenstates (ADS method) [6], and selfconjugate multibody states (BPGGSZ method) [7][8][9]. The most widely used D decays for the BPGGSZ method are D → K 0 S h + h − , where h = π, K. Measurements of γ using these final states have been performed by the Belle, BaBar and LHCb Collaborations [9][10][11]. Recently the first constraints on γ using the BPGGSZ method with a four-body D decay were reported [12]. BPGGSZ analyses require an understanding of the interference effects between D 0 andD 0 decays, especially concerning the strong-phase difference between the D 0 and D 0 decay amplitudes.
A precise measurement of the strong-phase difference in D → K 0 S,L π + π − decays was reported by the BESIII Collaboration recently [13]. The first measurements of the strong-phase difference between D 0 andD 0 decaying to the K 0 S,L K + K − final state were reported by the CLEO Collaboration, using a data set equivalent to an integrated luminosity of 818 pb −1 that was collected at a center-of-mass energy corresponding to the mass of the ψ(3770) resonance [14]. In this paper, we present an improved measurement of the strong-phase parameters for D → K 0 S,L K + K − decays, using a ψ(3770) data sample corresponding to an integrated luminosity of 2.93 fb −1 recorded by the BESIII detector. This measurement can be used as an input to the modelindependent measurement of γ using the BPGGSZ method. Moreover, these strong-phase parameters serve as an essential input to the model-independent determination of charmmixing parameters and in probing CP violation with D → K 0 S,L K + K − decays [15]. The D → K 0 S K + K − decay proceeds via various intermediate resonances, which leads to a significant strong-phase variation over the phase space. We define the kinematic variables m 2 ± = (P K 0 S + P K± ) 2 , which serve as the basis of the D → K 0 S K + K − Dalitz plot. Here, P i (i = K 0 S , K + , K − ) is the four-momentum of the D decay product. The amplitude for B − → D(K 0 S K + K − )K − at (m 2 + , m 2 − ) can be written as where r B is the ratio of the magnitude of the suppressed to the favored B-decay amplitude, δ B is the CP -conserving strongphase difference between favored and suppressed B-decay amplitudes, γ is the weak-phase difference between the B decay amplitudes, and f D (m 2 + , m 2 − ) fD(m 2 + , m 2 − ) is the am-plitude of the D 0 → K 0 S K + K − (D 0 → K 0 S K + K − ) decay. We neglect CP violation in D decays as in Ref. [8], and can thus use the relation fD(m 2 + , m 2 − ) = f D (m 2 − , m 2 + ) so that Eq. (1) can be written as Therefore, the decay rate of a B − meson is where ∆δ D ≡ δ D (m 2 + , m 2 − ) − δ D (m 2 − , m 2 + ), and δ D (m 2 + , m 2 − ) is the strong phase of f D (m 2 + , m 2 − ). Hence, knowledge of ∆δ D is essential for the determination of γ in B − → DK − decays.
In the literature, both model-dependent and modelindependent BPGGSZ methods are used. In the modeldependent approach, the D 0 amplitude is obtained using a flavor-tagged D 0 meson sample selected from the D * ± → D 0 (K 0 S K + K − )π ± decay, which is fit to an amplitude model describing the decay of D 0 → K 0 S K + K − [16] to determine f D (m 2 + , m 2 − ). The amplitude model is then used in an unbinned likelihood fit to the B-meson data sample to determine γ, δ B , and r B . However, this method results in a model-dependent systematic uncertainty on the measured value of γ which is difficult to quantify [17]. These modeldependent uncertainties have been estimated to lie between 3 • to 9 • [18,19], which limits the precision on γ that future measurements performed with much larger B-meson data samples [20,21] can obtain.
An alternative method of measuring γ is in a modelindependent manner that relies on defining a number of bins in the D 0 → K 0 S K + K − Dalitz plot [8]. This approach determines γ from the measured rate in each bin of the Dalitz plot, rather than fitting the Dalitz plot distribution to an amplitude model. The method requires information about ∆δ D (m 2 + , m 2 − ) in each bin, which is accessible at the ψ(3770) resonance by exploiting the quantum coherence of the D 0D0 pair produced in ψ(3770) decays. The advantage of this method is that the hard-to-quantify systematic uncertainty related to the model assumption is replaced by the uncertainty on the binned strong-phase parameters of the D decay mode. These strong-phase parameter uncertainties are statistically dominated, and thus well understood. The major disadvantage of the model-independent method is the inevitable loss of information that arises from binning, which reduces the statistical sensitivity of the γ measurement by approximately 80% compared to the model-dependent method [14].
The remainder of this paper is structured as follows. In Sec. II we define the formalism used to measure the strongphase parameters with ψ(3770) data. We explain the Dalitzplot binning in Sec. III. In Sec. IV we outline the features of the BESIII detector and the simulation techniques used in the analysis. We describe the event-selection criteria and the procedure for estimating the data yields in Secs. V and VI, respectively. In Sec. VII we explain the procedure for estimating the bin yields, including the various corrections applied. We describe the extraction of strong-phase parameters and the calculation of systematic uncertainties in Secs. VIII and IX, respectively. We present a discussion on the impact of these results on γ in Sec. X. In Sec. XI we give the conclusion and outlook.

II. FORMALISM
The model-independent method [8] for a three-body D decay is implemented as follows. The entire Dalitz plot is divided into 2N bins, with N bins symmetrically placed on either side of the m 2 + = m 2 − line. We follow the convention in which bins with m 2 + ≥ m 2 − are labelled with i and bins with m 2 + < m 2 − are labelled with −i. Thus, the 2N bins are assigned labels from −N to N excluding zero. The interchange of the Dalitz plot variables m 2 + ↔ m 2 − corresponds to the interchange of positions of the bins i ↔ −i. In order to extract the strong-phase difference parameters, we need to determine the yield in each bin for flavor-, CP -and mixed-CP tagged D → K 0 S K + K − decays. The number of flavor-tagged D 0 → K 0 S K + K − decays K i in the ith bin of the Dalitz plot is defined as where a D is a normalization factor equal to the total number of D 0 → K 0 S K + K − decays in the flavor-tagged charm sample, F i is the fraction of D 0 → K 0 S K + K − decays in the ith bin, and the integral is over the (m 2 + , m 2 − ) region defined by the ith bin. Here and elsewhere the values of K (−)i are corrected for efficiency and also for the presence of any doubly Cabibbo-suppressed (DCS) component (see Sec. VII A). We assume f D (m 2 + , m 2 − ) has been normalized such that where the integral is over the whole Dalitz plot. For each bin the interference between D 0 andD 0 decays can be parameterized by two variables c i and s i , which are the amplitude-weighted averages of cos ∆δ D and sin ∆δ D , defined as: and From Eqs. (6) and (7) it is evident that The equality is satisfied only if f D is constant throughout the bin. Thus the yield of B ± decays in the ith bin, N i , is obtained by integrating Eq. (3), which results in A maximum likelihood fit to binned B − → DK − decay yields, using Eq. (8) as a probability density function with externally measured values of c i and s i as inputs, then allows γ to be determined along with r B and δ B .
We now describe how ψ(3770) data are used to determine the values of c i and s i . The D 0D0 pair from the decay of the ψ(3770) (or if directly produced from the virtual photon in an e + e − annihilation) is in a C-odd eigenstate, as long as there are no additional particles in the final state. This quantum correlation between the mesons leads to the total D 0D0 decay rate being sensitive to the strong-phase difference between the D 0 andD 0 amplitudes. For example, the decay of one D to a CP -even eigenstate fixes the other D to the CP -odd admixture of D 0 −D 0 / √ 2. Hence, if the other D decays to K S,L K + K − , the total rate will be sensitive to the interference between the D 0 andD 0 amplitudes and the strong-phase parameters. Generally, this interference affects the decays of one D in combination with the other. If only one D meson is reconstructed, leaving the companion D meson to decay to any final state, the decay rate is largely insensitive to the effects of quantum correlations; we refer to the reconstructed samples of such events as single-tag (ST) decays. If both D mesons are required to be in specific final states, the rates can be significantly enhanced or suppressed in the quantum-correlated events compared to the expected rate if the decays are uncorrelated; we refer to the reconstructed samples of such events as double-tag (DT) decays. Hereafter, all the D decay final states, except the signal mode K 0 S,L K + K − , are referred to as "tags".
The D → K 0 S K + K − decay amplitude from a CP eigenstate is where + (−) indicates a CP -even (CP -odd) state. Therefore, the expected number of events M ± i in the ith bin of a sample that has been tagged with a decay that has a CP -even fraction F + is where S ± (S f ) are the efficiency-corrected single-tag yields of the CP -eigenstate (flavor) modes used in the analysis and DT,i is the DT efficiency in the ith bin. The value of F + is equal to 1 (0) for a pure CP -even (CP -odd) tag mode. We refer to modes with intermediate values of F + as quasi-CP tags. The values of c i alone can be extracted using Eq. (10). The relation s i = 1 − c 2 i is a good approximation only for N > 200 [22], which is not feasible with the available data sample. However, analysing D → K 0 S K + K − decays tagged by D → K 0 S h + h − (h = π, K) decays gives access to both c i and s i . The amplitude of the D 0D0 pair produced by the ψ(3770) decaying to K 0 where (m 2 + , m 2 − ) are the Dalitz plot coordinates corresponding to the phase space of the K 0 S h + h − decay. The expected event rate in which one D decays in the region of phase space defined by the ith bin of the D → K 0 S K + K − Dalitz plot and the other D in the region of phase space defined by the jth bin of the D → K 0 S h + h − Dalitz plot can be written as where N D 0D0 is the number of D 0D0 pairs in the ψ(3770) data sample and DT,ij is the DT efficiency in the ith and jth pair of bins. The two-fold ambiguity in the sign of s i can be resolved using weak amplitude-model assumptions. Note that Eq. (10) is symmetric under the interchange of i ↔ −i and Eq. (12) is symmetric under the interchange of pair, (i, j) ↔ (−i, −j) and (i, −j) ↔ (−i, j). No such symmetry exists for the values of K i because f D (m 2 + , m 2 − ) = f D (m 2 − , m 2 + ). Ignoring the very low level of CP violation in the neutral kaon system, the K 0 state is an equal admixture of K 0 S and K 0 L states. Therefore, in the decays of correlated D 0D0 pairs we expect a significant fraction of the D mesons to decay to the K 0 L K + K − final state as well. Although so far γ has only been determined using D → K 0 S h + h − decays, the decay D → K 0 L K + K − has a close connection with D → K 0 S K + K − that can be exploited to improve the precision with which c i and s i are determined. In the absence of CP violation, CP |K 0 S = |K 0 S and CP |K 0 Therefore, the number of events in the ith bin of a CP -or quasi-CP tagged D → K 0 where K i is defined in analogous fashion to K i (see Eq. (4)). Furthermore, the expected event rate in the ith bin of the D → K 0 L K + K − Dalitz plot and the jth bin of the D → K 0 S h + h − Dalitz plot can be written as The symmetries between the exchange of coordinates in the cases of M i and M ij are also present for M i and M ij .
In order to improve the precision of the extracted values of c i and s i constraints are imposed on the difference ∆c i = c i − c i and ∆s i = s i − s i ; these constraints are explained in Sec. VIII.
All the relations given in Sec. II are independent of the shape of the Dalitz plot bins. The original proposal [8] was to divide the Dalitz plot into rectilinear bins. The reduction in sensitivity of such an approach compared to an unbinned analysis is about 30% even with 20 bins [22]. The sensitivity of the model-independent method as a function of the bin shape is discussed in Ref. [22]; this paper concludes that binning schemes that minimise the variations of ∆δ D within each Dalitz plot bin give significantly improved statistical sensitivity compared to the rectilinear binning. An amplitude model can be used to guide the definition of bin boundaries in order to minimize the ∆δ D variation. The number of bins that can be used in the analysis is restricted by the available statistics in either the ψ(3770) or B-decay data samples. Since the use of an amplitude model is only to define the bin shapes, the model neither leads to any bias nor introduces any model-dependent uncertainties on the measurement of γ. However, a model that poorly describes the phase variation of the amplitude over the Dalitz plot may lead to a lower than expected statistical sensitivity to γ.
In the current analysis we employ an amplitude model for D 0 → K 0 S K + K − decays developed by the BaBar Collaboration [16] to define the bin shapes. Our choice of model and bin definitions is consistent with the previous measurement [14]. The amplitude model is constructed in the isobar formalism, where the amplitude at a phase-space point is defined as a coherent sum of two-body amplitudes and a non-  (16) as shown in Fig. 1 for N = 2, 3 and 4. Here i = 1, 2, ..., N are the bin numbers. The bins in the region m 2 + > m 2 − are defined symmetrically. The class of binning defined by Eq. (16) is referred to as the "equal-∆δ D " binning scheme. A smaller number of bins is the best choice to measure c i and s i precisely, but this will potentially reduce the sensitivity to γ. On the other hand, a larger number of bins provides increased sensitivity to γ, because it is a better approximation to the unbinned method. Keeping this trade-off in mind, we perform the analysis with N = 2, 3 and 4 bins. Binning the Dalitz plot with N > 4 is not yet feasible with the size of the sample of ψ(3770) data collected by BESIII; the fit to determine c i and s i described in Sec. VIII fails with N > 4.
In order to ascertain the quality of the binning, a figure-ofmerit based on the ratio of statistical sensitivity of the binned to the unbinned approach, known as the binning quality factor, Q, is defined in Ref. [22]. The predicted values of Q for this model are determined to be 0.771, 0.803 and 0.822 for N = 2, 3 and 4 bins, respectively [14]. The measured values were 0.94 +0. 16 −0.06 , 0.87 +0.14 −0.06 and 0.94 +0.21 −0.06 for N = 2, 3 and 4 bins respectively [24]. Since these values are close to one it implies that the loss of sensitivity due to the current bin definitions is small. An optimal binning scheme, which accounts for the distribution of the B-meson data sample across the Dalitz plot, as well as the ∆δ D variation, is found to give negligible improvement to the projected sensitivity compared to the "equal-∆δ D " binning [14]; hence, it is not pursued further.

IV. BESIII DETECTOR AND EVENT GENERATION
We analyse an e + e − collision data sample produced by the Beijing Electron Positron Collider II (BEPCII), which corresponds to an integrated luminosity of 2.93 fb −1 [25], collected by the BESIII detector at a center-of-mass energy of √ s = 3.773 GeV. The BESIII experiment is a general purpose solenoidal detector with a geometrical acceptance of 93% of the 4π solid angle. It has a He-gas-based multilayer drift chamber (MDC) for measuring the momentum and specific ionization loss (dE/dx) of the charged particles, a plasticscintillator-based time-of-flight (TOF) measurement system for the identification of charged particles, and an electromagnetic calorimeter (EMC) consisting of CsI(Tl) crystal, which is used to measure the energy of the neutral showers and identify electrons. The detector is encapsulated in a magnetic field of 1 T provided by a superconducting solenoid. A resistiveplate-chamber-based muon counter is interleaved between the flux-return yoke of the magnet. The MDC has a transversemomentum resolution of 0.5% at 1 GeV/c. The time resolution of the TOF is about 80 ps in the barrel region and 110 ps in the endcap region, enabling a 2σ K/π separation up to a momentum of 1 GeV/c. The energy resolution of the EMC for 1 GeV photons is about 2.5% in the barrel region and 5% in the endcap regions. More details about the BESIII detector can be found in Ref. [26].
Simulated samples produced with the GEANT4-based [27] Monte Carlo (MC) package, which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiency and to estimate the backgrounds. The simulation includes the beam-energy spread and initial-state radiation (ISR) in the e + e − annihilations modelled with the generator KKMC [28]. The inclusive MC samples consist of the production of DD pairs, the non-DD decays of the ψ(3770), the ISR production of the J/ψ and ψ(3686) states, and the continuum processes incorporated in KKMC [28]. The known decay modes are modelled with EVTGEN [29] using branching fractions taken from the Particle Data Group [30], and the remaining unknown decays from the charmonium states with LUNDCHARM [31]. The final state radiation (FSR) from charged final-state particles is incorporated with the PHOTOS package [32]. The simulation of quantum-correlations in the process ψ(3770) → D 0D0 is done outside the EVTGEN framework, using an algorithm developed by the CLEO Collaboration [33]. The effective integrated luminosity of the generated D 0D0 sample is about ten times that of the data. For the efficiency determination we use signal MC samples. Signal MC samples consist of D 0 → S tag ,D 0 → X decays for the reconstruction of STs and D 0 → K 0 S,L K + K − ,D 0 → S tag decays for the recon-struction of DTs, where S tag is a tag final state and X is any inclusive final state. Each signal MC sample corresponds to a specific ST or DT decay mode studied in this paper and contains 2 × 10 5 events.

V. EVENT SELECTION
In this section we initially describe the requirements for selecting the reconstructed particles that are combined to form the final states of interest. Then we present the selection criteria of fully reconstructed tag modes and partially reconstructed tag modes in Secs. V A and V B, respectively. Table I summarizes the set of tag modes used to reconstruct D 0 final states. The decay channels are split into five categories: Signal, Flavored, Mixed CP , CP -odd and CPeven. A highlight of this analysis is that the quasi-CP mode D → π + π − π 0 , which has a large branching fraction, is used for the first time for the strong-phase measurements in the D → K 0 S,L K + K − analysis. The F + value of π + π − π 0 is measured in Ref. [34] and the mode is found to be overwhelmingly CP -even. Hence in this analysis we treat π + π − π 0 as a CP -even tag taking into account its F + value. In the analysis, daughter particles are reconstructed as: In this section we will describe the selection criteria implemented to reconstruct these final states.

Type
Tag modes Signal For the charged tracks the polar angle θ is required to be within the MDC acceptance, which is |cos θ| < 0.93. The distance of closest approach of a primary track from the interaction region is required to be less than 10 cm in beam direction and less than 1 cm in the plane perpendicular to the beam direction to remove tracks not originating from e + e − collisions. For neutral showers the energy deposited in the EMC is required to be larger than 0.025 GeV in the barrel region (|cos θ| < 0.8) and larger than 0.050 GeV in the endcap region (0.86 < |cos θ| < 0.92), which reduces the effect of electronic noise and deposits resulting from beam-related backgrounds. Moreover, the angle between the position of the shower and any extrapolated charged track in the EMC must be greater than 10 • to reduce the number of showers related to charged tracks. Furthermore, we require the time of the shower to be less than 700 ns after the event start-time to further suppress fake photons associated with electronic noise and beam backgrounds.
The particle identification (PID) is performed by combining the dE/dx information from the MDC as well as the time-offlight of the charged particle. The likelihoods for the kaon hypothesis L K and pion hypothesis L π are calculated. Tracks satisfying the condition L K > L π are identified as kaons and vice versa for pions. For electrons the PID is performed by defining a likelihood based on information about dE/dx in the MDC, time-of-flight and deposited energy and shape of the electromagnetic shower from the EMC. The track is identified as an electron if L e /(L e + L K + L π ) > 0.8 and L e > 0.001, where L e is the likelihood of the electron hypothesis.
A K 0 S candidate is formed by considering a pair of intersecting oppositely charged tracks. These tracks are not subject to any track quality requirement or PID. The closest approach of these tracks to the interaction point is required to be less than 20 cm along the beam direction with no requirement in the transverse direction. A secondary vertex fit is performed to form the K 0 S vertex, and candidates with χ 2 < 100 are selected. The updated four momenta after the secondary vertex fit are used later in this analysis. The mass of a K 0 S candidate is required to be within the range (0.487, 0.511) GeV/c 2 . In order to suppress combinatorial backgrounds from two pions that are not from a true K 0 S , the flight significance, L/σ L , is required to be greater than two, where L is the flight length and σ L is the uncertainty in L from the secondary vertex fit.
Both π 0 and η candidates are reconstructed from a pair of photons, where at least one of the photons must be reconstructed in the barrel region; this requirement reduces combinatorial backgrounds that arise from the large number of showers in the endcap region that are related to beam backgrounds. The invariant mass of the two photon candidates must be in the range (0.110, 0.155) GeV/c 2 or (0.480, 0.580) GeV/c 2 for π 0 and η candidates, respectively. In order to improve the momentum resolution, a kinematic fit of the two photons is performed with their invariant mass constrained to the nominal mass of π 0 or η meson taken from the PDG [30]. Only π 0 and η candidates with χ 2 < 20 are selected. The improved values of the momenta are used later in the analysis. For ω candidates the invariant mass of the π + π − π 0 combination is required to be within the range (0.760, 0.805) GeV/c 2 and for η candidates the invariant mass of the π + π − η combination is required to be within the range (0.938, 0.978) GeV/c 2 . All the invariant mass intervals described correspond to approximately ±3 times the standard deviation about the mean of the reconstructed distribution.

A. Selection of fully reconstructed tags
Fully reconstructed tags are decay modes that do not contain an undetected particle in the final state. Before describing the kinematic variables used to select fully reconstructed tags, we introduce two additional vetoes that remove specific backgrounds to certain tag modes. The first veto is to suppress backgrounds arising from cosmic rays and lepton-pair events in the ST reconstruction of the two-body decay channels K + K − , π + π − and K ± π ∓ . Here, we reject events in which the two charged tracks that reconstruct the ST candidate are consistent with being an e + e − or µ + µ − pair. In addition, to suppress cosmic muons, we reject events where the time-offlight difference between the two tracks is greater than 5 ns. Further, an event that has neither an EMC shower with an energy greater than 50 MeV nor an additional charged track in the MDC is rejected. The second veto is to remove the CPodd K 0 S π 0 , K 0 S → π + π − background to the predominantly CP -even π + π − π 0 tag mode; here we reject events that satis- S refers to the nominal mass of the K 0 S meson given in Ref. [30].
For all fully reconstructed tag modes, the selected finalstate particles are combined to reconstruct the D decay. Since the DD pair production occurs at the ψ(3770) resonance, there are no additional particles in the final state, so the energy of each D meson is equal to √ s/2. Thus, with a well measured beam energy E beam (= √ s/2) we define two quantities to reconstruct the D candidates: the energy difference, and the beam-constrained mass, Here E i and p i are the energies and momenta of the D decay products in the center-of-mass frame. Properly reconstructed candidates will peak at zero in the ∆E distribution and at the nominal mass of the D 0 meson [30]  Partially reconstructed tags collectively refer to the tag modes where there is one particle in the final state, either a K 0 L meson or a neutrino, which is not reconstructed. Modes with more than one missing particle in the final state are not considered in this analysis. Due to the presence of a missing particle in the final state, these tag modes can be reconstructed only as DTs so that four-momentum conservation can be exploited in the reconstruction.
Selections of partially reconstructed tag modes proceed as follows. The opposite-side D candidate is reconstructed as a ST candidate using the criteria given in the Sec. V A. All the particles except the missing particle in the final state are reconstructed from the unused tracks and showers that satisfy the selection criteria already described. The presence of an unreconstructed K 0 L is inferred from the missing-mass distribution, calculated from the missing energy, E miss , and the missing momentum, p miss , in the center-of-mass frame as which peaks at m 2 K 0 for signal, where m K 0 is the nominal mass of the neutral kaon given in Ref. [30]. The presence of a neutrino is inferred using the variable which peaks at zero for signal. Again we take advantage of resonant production and the knowledge of beam energy to determine E miss and p miss . Figure 2 shows example distributions of M 2 miss and U miss . Reconstruction using the missingmass technique inevitably results in a higher level of background than the full-reconstruction method. To reduce the background further, we do not consider events that have more charged tracks or neutral particles than required in the final state. The angle, α, between the p miss and the nearest unassigned shower is calculated. All the events with cos α > 0.98 are retained. For the events with cos α < 0.98 modedependent criteria are applied on the energy of the unassigned shower. Even though we reject events with additional neutral particles in the final state, there is significant background in the modes with neutral particles in the final states, arising from the final states having additional neutral particles that are not reconstructed. For example, in the case of D → K 0 L π 0 decays there are backgrounds from K 0 L π 0 π 0 where one π 0 meson is not reconstructed, so the event passes all our selection criteria. These backgrounds can be further reduced by applying criteria on the momentum spectrum of reconstructed π 0 or η candidates wherever applicable. The values of these criteria are selected based on optimization studies that use the inclusive MC samples. This optimization maximizes the figure-ofmerit defined as S/ √ S + B, where S (B) are the number of signal (background) events in the signal region retained by the selection; the signal region for the optimization is the interval 0.2 < M 2 miss < 0.3 GeV 2 /c 4 . The values of the shower energy and π 0 (η) momentum criteria are varied, and the value that maximizes the figure-of-merit is chosen.

VI. ESTIMATION OF ST AND DT YIELDS
In Secs. VI A and VI B we will describe the methods of estimating ST and DT yields, respectively. Note that DT yields are only required bin-by-bin, not integrated over the Dalitz plot.

A. ST yields
ST yields of fully reconstructed tag modes are determined from maximum likelihood fits to the M BC distribution. Our probability density function (PDF) is a sum of the signal shape derived from the signal MC sample convolved with a Gaussian function to account for any difference in resolution between data and MC simulation, and an ARGUS function [35] to model the background. The threshold of the ARGUS function is fixed at M BC = 1.8865 GeV/c 2 , which corresponds to the kinematic limit of D 0 production at the ψ(3770). The peaking background is modelled by the shapes and yields obtained from the inclusive MC sample; this assumption is considered as a source of systematic uncertainty. The flavor-tag modes D 0 → K − π + and D 0 → K − π + π 0 have a peaking background of approximately 0.2% from DCS decays. The dominant peaking background to the decays D → K 0 S π 0 and D → K 0 S π 0 π 0 is from D → π + π − π 0 (0.5%) and D → π + π − π 0 π 0 (7%) decays, respectively. The M BC distribution is fit over the range (1.83, 1.88) GeV/c 2 . The ST yields are obtained by integrating the M BC distribution in the range (1.86, 1.87) GeV/c 2 . In order to eliminate the small effect of D 0D0 mixing, the measured ST yields of CP modes are multiplied by a correction factor of 1/(1 − η ± y D ), where η ± is the CP eigenvalue of the mode and y D is the charm-mixing parameter taken from Ref. [30].
The ST yield, S ST , of a partially reconstructed tag is calculated using the relation where N D 0D0 is the number of D 0D0 pairs in the BESIII data sample [36] and B ST is the branching fraction of the tag mode, which is taken from Ref. [30] where available. The branching fractions of all D → K 0 L X modes except D → K 0 L π 0 are not available in Ref. [30], hence we assume the branching fractions of these modes to be the same as for the corresponding D → K 0 S X modes. We note that this reasoning is not strictly valid, as the interference between Cabibbo-favored (CF) (D 0 →K 0 X) and DCS transitions (D 0 → K 0 X) can lead to a difference in the decay rates for D → K 0 L X and D → K 0 S X. However, this difference is expected to be less than 10% [37], which is considered as a systematic uncertainty; the difference will barely affect our final results, as the ST yields are used only for yield normalisation, as given in Eqs. (10) and (14). The ST yields calculated using Eq. (21) have larger uncertainties compared to the fully reconstructed tags, largely due to the uncertainty of the assumed values of B ST . The ST M BC fits are shown in Fig. 3 and the yields are given in Table II. The effect on the final measurement due to the uncertainty in the measured values of the ST yields is treated as a systematic uncertainty. The ST yield uncertainty includes systematic uncertainties related to the fit procedure.

B. DT yield
For fully reconstructed DT modes we follow a sidebandestimation method developed by the CLEO Collaboration [38] to determine the DT yield. The sidebands are defined on the two-dimensional (M D 0 BC , MD 0 BC ) plane as shown in Fig. 4. Here, the M D 0 BC (MD 0 BC ) refers to the M BC distribution of signal (tag) side. In Fig. 4, S refers to the signal region, sideband A (B) contains events which are from misreconstructed tag (signal) decays, sideband C consists of continuum events and sideband D consists of events that are purely combinatoric. The amount of combinatorial (non-peaking) backgrounds in the signal region is estimated from the events in the sideband regions. Thus the total DT yield, N DT , of K 0 S K + K − is estimated as where a i is the area of the corresponding region i = A, B, C, D or S, N i refers to the yields in the sideband region, N S is the yield in the signal region before background correction (uncorrected yield) and N P is the peaking-background yield estimated from the MC simulation (see Sec. VII B).
In the case of partially reconstructed tag modes we follow a similar sideband-estimation method as in Ref. [24]. Here three regions are defined on the M 2 miss or U miss distributions: low sideband (L), signal region (S) and high sideband (H). The total yield is estimated as where N S , N L and N H are the uncorrected yields in the signal and sideband regions, N P i refers to the peaking background in the ith region, δ and γ refer to the ratio of combinatorial backgrounds in the signal region to that in the L and H sideband regions, respectively, and α and β refer to the ratio of signal in region S to that in the regions L and H, respectively. The values of α, β, δ and γ are derived from MC samples. Here the definitions of sidebands are mode dependent. We follow the same optimization procedure described in Sec. V B to define the signal regions. The peaking backgrounds are estimated from MC samples as described in Sec. VII B. pure CP eigenstates and mixed CP states; we highlight the important differences.
In order to improve the resolution on the Dalitz plot variables (m 2 + , m 2 − ), a kinematic fit is performed for the two pions from the K 0 S candidate obtained after the secondary vertex fit are combined with the h + and h − into a common fit to the nominal mass of the D 0 meson taken from Ref. [30]. In the case of a D → K 0 L h + h − candidate, a missing particle is created using the position of an EMC shower associated with the K 0 L candidate. The mass of this object is set to the nominal mass of the K 0 L meson taken from Ref. [30]; it is combined with h + h − tracks and fit to the nominal mass of the D 0 meson. A 35 to 40% (30 to 35%) improvement in the m 2 ± reso-   The Dalitz plot distribution of the D → K 0 S K + K − candidates reconstructed against CP -even tag modes and their corresponding M 2 K + K − projections are given in Fig. 5. The presence of a significant peak around M 2 K + K − ∼ 1.04 GeV 2 /c 4 is due to the decay D 0 → K 0 S φ, φ → K + K − . These events are distributed along the diagonal boundary of the Dalitz plot. As D 0 → K 0 S φ constitutes a large fraction of the total D 0 → K 0 S K + K − decay width [30], a higher population of events is seen in the region enclosing the φ resonance. A similar peak is absent in the M 2 K + K − distribution of D → K 0 S K + K − candidates reconstructed against CPodd tag modes shown in Figs. 5(c) and (d). This is a consequence of the quantum-correlation in data. Since each D meson is of opposite CP eigenvalue, the K 0 S K + K − candidates reconstructed against CP -odd tags should decay through CPeven intermediate states. Hence it cannot decay through the D → K 0 S φ state. The dominant CP -even intermediate state is the D → K 0 S a(980) 0 decay. The distribution of events in the Dalitz plot is observed to be flatter than in the case of K 0 S K + K − tagged against a CP -even state. Since K 0 S and K 0 L have opposite CP eigenvalues, the entire scenario is reversed in the case of D → K 0 L K + K − decays as shown in Fig. 6. The Dalitz plot distribution of D → K 0 L K + K − candidates against CP -even modes resembles that of D → K 0 S K + K − candidates against CP -odd modes and vice versa. against the self-conjugate mode D → K 0 S π + π − is given separately for signal and tag sides in Fig. 7. The Dalitz plot of D → K 0 S π + π − tags is consistent with that presented in Ref. [13]. In Fig. 7(d), the enhancement above M 2 π + π − ∼ 1.3 GeV 2 /c 4 corresponds to D → K * (892) ± π ∓ decays, whereas the peak around M 2 π + π − ∼ 0.6 GeV 2 /c 4 corresponds to D → K 0 S ρ 0 decays. The D → K * (892) ± π ∓ decays can be seen as two bands that are parallel to the vertical and horizontal axes of the Dalitz plane. The decay D → K 0 S ρ 0 lies close and parallel to the diagonal boundary. Since the D → K 0 L K + K − decays reconstructed against D → K 0 S π + π − decays are not in a CP eigenstate, the Dalitz plot distribution is a combination of both the CP -even and CP -odd tagged K 0 L K + K − Dalitz plots. The Dalitz plot structure of D → K 0 S K + K − reconstructed against D → K 0 S,L π + π − has similar features to those shown in Fig. 7.

A. Dalitz plot binning, bin yield estimation and corrections
In this section we describe our method of binning the Dalitz plots and calculating the bin yields and efficiencies. The pro-cedures for correcting the bin migration and DCS correction for flavor-tag yields are also explained.
The binning prescription followed in our analysis is described in Sec. II. The entire D 0 → K 0 S,L K + K − Dalitz plot is divided into N = 2, N = 3 and N = 4 equal-∆δ D bins. In the case of D → K 0 S,L π + π − tag modes, the entire Dalitz plot is divided into N = 8 equal-∆δ D bins identical to those defined in Ref. [13]. The uncorrected bin yields are obtained by counting the number of events in each bin. The bin yield M ± i (see Eq. (10)) of D → K 0 S K + K − reconstructed against CP tags and the flavor-tag yield, K i (see Eq. (4)) are calculated separately for each mode. The events in the ith bin of the D → K 0 S K + K − Dalitz plot and the jth bin of the D → K 0 S h + h − Dalitz plot are counted to obtain M ij (see Eq. (12)). A similar procedure is followed to obtain the yields K i , M ± i and M ij (see Eqs. (14) and (15)) for the D → K 0 L K + K − decay. The flavor-tag yield for the D 0 → K 0 S,L π + π − mode is taken from Ref. [13]. The yields of D → K 0 S K + K − decays reconstructed against CP tags are quite low. The inclusion of the D → π + π − π 0 tag mode results in an approximately 50% increase in the CP -even tag yield. The uncorrected yields of D → K 0 constructed against CP tags, along with their efficiencies, are given in Table II. Due to the finite (m 2 + , m 2 − ) resolution, events migrate between bins. Often these migrations are asymmetric between the bins because of the differing event densities in each bin. We correct for this using an unfolding method based on correction factors derived from the signal MC samples. For D → K 0 S,L K + K − decays reconstructed against CP and flavor tags, we define a 2N × 2N migration matrix U as where m ji are the events generated in the jth bin and reconstructed in the ith bin. The vector of migration-corrected data yields N and the vector of reconstructed yields in the signal region N S are related by In the case of D → K 0 S,L K + K − reconstructed against the D → K 0 S,L h + h − tags, the correlation between the bins on the signal and tag sides needs to be taken into account. Hence the total migration matrix is a tensor (Kronecker) product of signal-and tag-migration matrices. For a given number of signal bins N , the dimension of the migration matrix for The uncertainties in the matrix elements due to the finite size of the signal MC sample are considered as a source of systematic uncertainty. An example of the migration matrix for D → K 0 S K + K − candidates reconstructed against the D → K + K − tag mode is given in Table III. Typically the rate of migration out of bin 1, which contains the narrow φ resonance, is about 3% for D → K 0 S K + K − decays and about 5% for D → K 0 L K + K − decays. The rate of migration into bin 1 is significantly smaller due to the broader structures that occupy the remainder of the Dalitz plot away from the φ resonance. Throughout this unfolding procedure we assume signal and background migrate in identical fashion, because the background is dominated by peaking components.
The bin efficiency for each tag mode is evaluated from the signal MC sample. The signal MC yield in each bin is corrected for migration before calculating the efficiency. The bin efficiency is defined as the ratio of events reconstructed in each bin to the number of events generated. The bins are combined appropriately taking into account their symmetry (see Sec. II) when estimating the efficiencies. The total DT efficiencies are given in Table II. In the case of D → K 0 S K + K − , the efficiencies vary between (12.43 ± 0.07)% for K 0 π + π − distributions for D → K 0 S π + π − decay in the same events. whereas for D → K 0 L K + K − the efficiency varies between (15.85 ± 0.08)% for K 0 L K + K − vs. K − π + tags and (3.40 ± 0.04)% for K 0 S K + K − vs. K 0 L K + K − tags. The uncertainty on the efficiency is related to the size of the MC sample. The bin efficiencies are used to calculate the expected yield for each tag mode as given in Eqs. (10), (12), (14) and (15).
Both pseudo-flavor DT yields with F ∈ (K − π + , K − π + π 0 ) have contamination from DCS decays whose contribution is enhanced compared to ST yields due to the quantum correlation between the D 0D0 . Since these decays are used to determine K ( ) i , the presence of this DCS contamination may bias the values [39]. In order to correct for this effect, the yield in each bin is multiplied by a correction factor estimated using the decay model reported in Ref. [16]. The correction factors where r F D is the ratio of the moduli of the DCS to CF am-plitudes, for example |A(  [42] 0.81 ± 0.06 [42] for K ± π ∓ , and δ D F is the strong-phase difference between the DCS and CF amplitudes. The coherence factor, R F for flavormode F , accounts for the dilution in interference effects that arises when integrating over the phase space of multi-body decays [40]. The values of the parameters used to determine the correction factors are listed in Table IV. The fraction of events in each bin F ( ) i , defined in Eq. (4), is given in Table V. The [24]. Good agreement with the predicted values [14] is observed for the results given in Table V. The uncertainties in the final result due to the correction factors are small and are treated as a systematic uncertainty. The DCS correction is not required for the D 0 → K − e + ν e flavor-tag.

B. Bin-by-bin background estimation
In this section we explain the method used to estimate the peaking background. The amount of combinatorial background in each bin is estimated from the sideband-estimation methods described in Sec. VI B.
The peaking backgrounds are identified from the inclusive MC samples using the tool described in Ref. [43]. The backgrounds to fully reconstructed tags are found to be negligible. However, all the D → K 0 L X modes contain backgrounds from D → K 0 S X modes, where the π 0 mesons from K 0 S → π 0 π 0 decays are not reconstructed, so that the K 0 S is treated as a missing particle. The D → K 0 L X and D → K 0 S X decays are of opposite CP , hence the distribution of background events across the Dalitz plot is not the same as that for signal events. The level of these backgrounds varies between 2 to 4% depending on the tag mode. The bin-by-bin background estimation using the inclusive MC sample is not reliable for two reasons. First, there can be a difference between the branching fraction in data and that assumed in the MC simulation. Second, the MC samples are not tuned to reflect the distributions of events across the Dalitz plot. Both these issues will result in an incorrect estimation of the binby-bin background. Hence, we use a combination of data and background MC samples to estimate the backgrounds.
Our method of peaking background estimation is as follows. We generate dedicated background MC samples corresponding to each type of background decay. The background retention efficiencies for these backgrounds are calculated for each bin. The expected yields are calculated using the values of c i and s i obtained from the previous measurement [14] through the relations given in Eqs. (10), (12), (14) and (15). The yields are multiplied by the retention efficiencies to obtain the expected background yields in data. Though most of the combinatorial background beneath the D → K 0 L X signal decays are from non-resonant D → K + K − π 0 π 0 decays, these decays only contribute approximately 2% of the background in the signal region. The contributions from continuum backgrounds and low-lying charmonium resonance decays are found to be negligible for all DT modes, which is expected given the tight kinematic criteria that can be imposed close to open-charm threshold. The expected background yields are not subtracted from signal yield but added to the expected signal yield in the fit as explained in Sec. VIII.

VIII. EXTRACTION OF ci AND si
The uncorrected yields in related bins are combined according to the symmetry relations described in Sec. II. This procedure reduces the number of degrees of freedom in the fit.
are taken from Ref. [13]. The values of c The values are those reported in Ref. [24].
Here, N is the expected migration-corrected yield in a particular bin whose measured yield is N . P (N, N ) is the Poisson probability of observing a yield N given the mean N : If M represents the expected signal yield and B represents the expected background then N = M + B . It is to be noted that in Eq. (28) the comparison is between the uncorrected yield and the sum of expected signal and expected background in each bin. This method eliminates the possibility of unphysical negative bin yields arising from the subtraction of backgrounds from bins having low yields. The χ 2 term in Eq. (28) constrains the difference between the extracted values of c i (s i ) and c i (s i ) to lie within the uncertainties of the predicted differences ∆c i (∆s i ). The χ 2 is defined as is the predicted difference from the BaBar model [16] and σ ∆ci (σ ∆si ) is the uncertainty on ∆c i (∆s i ). The values of ∆c i (∆s i ) and σ ∆ci (σ ∆si ) are given in Table VI.
Large values of ∆c i and ∆s i are observed in bins i = 3 and i = 4 in the N = 3 and N = 4 binnings, respectively. These bins correspond to the lobes on the Dalitz plot that encompass the neutral resonance a 0 (1450) 0 . The model has very small amplitudes in this region, hence a small difference between the D → K 0 S K + K − and D → K 0 L K + K − models is proportionally more significant. Consequently, the uncertainties are also large. Improvement in the precision due to the χ 2 term is more significant for s i than c i . The minimization of Eq. (28) is performed using the MINUIT [44] package. The value of parameters obtained from the fit are given in Table VII. The fitting procedure is validated using pseudodata samples generated by a standalone simulation. The bin yields and backgrounds are generated according to the relations given in Sec. II. The effects of bin migration are also considered in the simulation. The bin yields are fluctuated assuming a Poisson distribution. The procedure is repeated 500 times to obtain a pull distributions for c The systematic uncertainties related to the ST yields are evaluated as follows. First, a combined systematic uncertainty in the yields due to various assumptions made in the fit is estimated. For example, the endpoints of the ARGUS function are fixed in our fits. We rerun the fits by changing the endpoints by ±0.5 MeV/c 2 , which is the approximate uncertainty related to the beam-energy spread in the endpoint. The difference between the new value of the yield and the nominal value is taken as a systematic uncertainty in the ST yield. Other assumptions include the estimation of peaking backgrounds from MC simulations and the assumed branching fractions in partially reconstructed modes. The statistical uncertainties in the yields from the fits are added in quadrature with the systematic uncertainties to obtain the total uncertainty related to the ST fit. The ST yield is smeared within the total uncertainty and the distribution of the resulting c ( ) i and s ( ) i values is obtained as laid out above. The systematic uncertainties due to the ST yield are found to be small. The systematic uncertainty due to the factors used to correct for the effect of charm mixing is found to be negligible. Furthermore, uncertainties related to the absolute efficiency do not affect the results due to the use of yield ratios and normalization constants fit to data.
The systematic uncertainty due to the finite statistics of the flavor-tag yields K i are evaluated by smearing the input values assuming a Gaussian distribution around the nominal value with width equal to the uncertainty. The systematic uncertainty due to the flavor-tag yield of the D → K 0 S π + π − decay is estimated using the covariance matrix reported in the Ref. [13]. The uncertainties due to flavor-tag yields are small due to the large yields compared to those for CP and D → K 0 h + h − tags.
In our fit the values of c ( ) i and s ( ) i of D → K 0 S,L π + π − are fixed to the values reported in Ref. [13]. The covariance matrix used to smear the values is also taken from Ref. [13]. This uncertainty is the dominant systematic uncertainty.
The value of the total number of D 0D0 pairs is fixed in our fits. This information is used to normalize the bin yields of D → K 0 S,L K + K − tagged by D → K 0 S,L h + h − decays. The related systematic uncertainty is evaluated by smearing the value of N D 0D0 within its uncertainty assuming a Gaussian distribution. Since the value of N D 0D0 pairs is measured precisely [36], the systematic uncertainty due to this input is negligible.
Signal MC samples are widely used in this analysis for various purposes such as constructing migration matrices and determining the selection efficiencies. The systematic uncertainties due to the finite size of the MC samples are evaluated by smearing the matrix elements within their uncertainties assuming a Gaussian distribution taking correlations into account. Any systematic uncertainty due to the incorrect MC modelling is cancelled since we use ratios of ST and DT yields in computing the values of M .
The effect of bin-by-bin peaking background from D → K 0 S h + h − decays in the D → K 0 L h + h − signal sample is estimated using the c i and s i values reported in Ref. [14]. To estimate the systematic uncertainty from the background parametrization, the estimated amount of background is varied by a Gaussian function within its statistical uncertainty. The level of the background with respect to the partially reconstructed tags is larger than that with respect to the fully reconstructed tags, hence the systematic uncertainties due to the background parametrization for partially reconstructed tags are larger than those for the fully reconstructed tags. Since the backgrounds are only identified but not estimated from MC simulations, there is no systematic uncertainty arising due to any difference in branching fractions between data and MC simulations.   [41]). In order to reduce the uncertainty due to the statistics of the B decay sample to a negligible level, the pseudodata sample size is set to 5×10 6 events. The bin yields are fluctuated according to a Poisson distribution to produce a new bin yield, N I . The expected bin yield, N E , is calculated using c i and s i values obtained by smearing the measured values by their uncertainties, assuming a Gaussian distribution and taking the correlations into account. The best fit values of r B , δ B and γ are extracted by minimizing The pseudodata simulations are repeated 10000 times and the resulting γ distributions for N = 2, 3 and 4 bins are shown in Fig. 9.
The distributions of all parameters (r B , δ B and γ) are found to be asymmetric. Using the root mean square (RMS) of the distribution, we estimate the uncertainties on γ to be 2.3 • , 1.3 • and 1.3 • for the schemes with N = 2, 3 and 4 bins, respectively. We also estimate an asymmetric uncertainty by integrating 16% of the distribution in the lower and upper tails to work out a 68% confidence level; the asymmetric uncertainties on γ are +2.4 • −1.8 • , +1.0 • −0.9 • and +0.9 • −0.9 • for schemes with N = 2, 3 and 4 bins, respectively. Better sensitivity for N = 3 and 4 compared to N = 2 is observed. This is due to the fact that dividing the data into more bins is a better approximation of the unbinned case. The lack of improvement in sensitivity while going from N = 3 to N = 4 is due to the nature of the D 0 → K 0 S K + K − Dalitz plot. The dominant resonances contributing to D 0 → K 0 S K + K − are D 0 → K 0 S φ and D 0 → K 0 S a(980) 0 , which are both located close to the M 2 K + K − kinematic limit. In a binning based on equal-∆δ D regions these bins are always enclosed by a similar pair of bins and the new pair of bins always encloses a region with a low density of events, as can be seen in Fig. 1.
Biases of up to 0.7 • are observed in the central values of γ for all binning schemes. A bias was reported in the previous analysis as well [14]. The bias is smaller with our measurement, but non-negligible given the size of the pseudodata sample used. In order to investigate the source of the bias, we rerun the pseudodata experiments ignoring any pairs of c i and s i values that lie outside the physical region given by c 2 i + s 2 i < 1; the bias is still observed, which is due to the non-Gaussian nature of the truncated c i and s i distributions. Therefore, instead of removing the unphysical values the uncertainties on the c i and s i are artificially reduced by a factor of three; in this case no observable bias is seen in any of the extracted parameters. Hence we conclude that the bias stems from some pairs of c i and s i values lying outside the physical region. Therefore, future measurements with a larger ψ(3770) sample [45] are likely to reduce or remove this bias.

XI. SUMMARY
Using a sample of ψ(3770) data corresponding to an integrated luminosity of 2.93 fb −1 collected by the BESIII detector, we report a measurement of the strong-phase difference parameters for D → K 0 S,L K + K − decays with the best precision to date. The results presented here are an important input to model-independent measurements of the CKM angle γ using the BPGGSZ method. These values are also essential for the model-independent determination of charmmixing parameters and in the search for CP violation in D 0 → K 0 S K + K − decays [46]. We note that the statistical uncertainties limit the precision of our measurements. Therefore, it is desirable to collect larger ψ(3770) data sets in future [45].
The major source of systematic uncertainty is due to the input strong-phase difference parameters of D → K 0 S,L π + π − decays [13]. Significant systematic uncertainties also arise from the background parametrization of the partially reconstructed D → K 0 L X decay modes. Both of these uncertainties depend on the size of the charm sample available.
Good agreement with the previous measurements by the CLEO Collaboration [14] is obtained in all bins. Hence, we have performed a combination of the BESIII and CLEO results, which is reported in App. C. The predictions from the BaBar model [16] lie within one to two standard deviations from our values. We have recently reported an amplitude model and branching fraction for D 0 → K 0 S K + K − decays