Measurement of the $D \to K^-\pi^+\pi^+\pi^-$ and $D \to K^-\pi^+\pi^0$ coherence factors and average strong-phase differences in quantum-correlated ${D\bar{D}}$ decays

The decays $D\to K^-\pi^+\pi^+\pi^-$ and $D \to K^-\pi^+\pi^0$ are studied in a sample of quantum-correlated $D\bar{D}$ pairs produced through the process $e^+e^- \to \psi(3770) \to D\bar{D}$, exploiting a data set collected by the BESIII experiment that corresponds to an integrated luminosity of 2.93 fb$^{-1}$. Here $D$ indicates a quantum superposition of a $D^0$ and a $\bar{D}^0$ meson. By reconstructing one neutral charm meson in a signal decay, and the other in the same or a different final state, observables are measured that contain information on the coherence factors and average strong-phase differences of each of the signal modes. These parameters are critical inputs in the measurement of the angle $\gamma$ of the Unitarity Triangle in $B^- \to DK^-$ decays at the LHCb and Belle II experiments. The coherence factors are determined to be $R_{K3\pi}=0.52^{+0.12}_{-0.10}$ and $R_{K\pi\pi^0}=0.78 \pm 0.04$, with values for the average strong-phase differences that are $\delta_D^{K3\pi}=\left(167^{+31}_{-19}\right)^\circ$ and $\delta_D^{K\pi\pi^0}=\left(196^{+14}_{-15}\right)^\circ$, where the uncertainties include both statistical and systematic contributions. The analysis is re-performed in four bins of the phase-space of the $D \to K^-\pi^+\pi^+\pi^-$ to yield results that will allow for a more sensitive measurement of $\gamma$ with this mode, to which the BESIII inputs will contribute an uncertainty of around 6$^\circ$.

Abstract: The decays D → K − π + π + π − and D → K − π + π 0 are studied in a sample of quantum-correlated DD pairs produced through the process e + e − → ψ(3770) → DD, exploiting a data set collected by the BESIII experiment that corresponds to an integrated luminosity of 2.93 fb −1 . Here D indicates a quantum superposition of a D 0 and aD 0 meson. By reconstructing one neutral charm meson in a signal decay, and the other in the same or a different final state, observables are measured that contain information on the coherence factors and average strong-phase differences of each of the signal modes. These parameters are critical inputs in the measurement of the angle γ of the Unitarity Triangle in B − → DK − decays at the LHCb and Belle II experiments. The coherence factors are determined to be R K3π = 0.52 +0. 12 −0.10 and R Kππ 0 = 0.78 ± 0.04, with values for the average strong-phase differences that are δ K3π

Introduction
In the Standard Model, CP violation is described by the irreducible complex phase of the Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix [1,2]. This matrix can be represented geometrically by the Unitarity Triangle in the complex plane, with angles α, β and γ (also denoted φ 2 , φ 1 and φ 3 ). The angle γ, equal to arg(−V us V ub * /V cs V cb * ) at O(λ 4 ) ∼ 10 −3 , is of particular interest, as it can be measured with negligible theoretical uncertainty in decays of the class B − → DK − , where D represents a superposition of D 0 andD 0 mesons reconstructed in a final state common to both. The γ sensitivity arises from the phase difference between the b → u and b → c amplitudes, which manifests itself in an interference term in the partial width. These decays are mediated by tree-level amplitudes and thus their rates are not expected to deviate from the Standard Model prediction.
In contrast, other observables that can be used to infer γ involve loop-level amplitudes, which are sensitive to non-Standard-Model processes. Comparison between the two sets of measurements therefore provides a powerful method to probe for new physics beyond the Standard Model. An important category of charm decays with which to perform B − → DK − measurements is D → K − nπ, where nπ signifies n ≥ 1 pions [3,4] 1 . Here the final state can be accessed either by the Cabibbo-favoured (CF) amplitudes of a D 0 decay, or the doubly Cabibbo-suppressed (DCS) amplitudes of aD 0 decay, or either of these amplitudes in conjunction with those responsible for D 0 -D 0 oscillations. In the case of final states where n ≥ 2, such as D → K − π + π + π − or D → K − π + π 0 , the interference term involving γ is also dependent on the variation of the D 0 andD 0 decay amplitudes over the multi-body phase space. Inclusive measurements that integrate over this phase space can be interpreted in terms of γ provided that three charm hadronic parameters are known. Generalising for the decay D → S, where S denotes a K − nπ final state, these parameters are the coherence factor R S , and the amplitude ratio r S D and CP -conserving strong-phase difference δ S D between the CF and DCS amplitudes averaged over phase space [5]: Here A S (x) is the decay amplitude of D 0 → S at a point in multi-body phase space described by parameters x and A 2 S = |A S (x)| 2 dx, with equivalent definitions for the D 0 →S amplitude. (The amplitudes are normalised such that, in the absence of oscillations, A 2 S corresponds to the branching ratio of the D 0 meson into mode S.) Maximum (zero) interference in the term sensitive to γ occurs in the limiting case R S = 1 (0). The parameter r S D is of the order λ 2 ≈ 0.06 for D → K − nπ decays.
Strong constraints exist on the values of r S D from measurements of branching fractions. The coherence factors and strong-phase differences may be determined through decays of quantum-correlated DD pairs produced in e + e − collisions at the ψ(3770) resonance. The decay rate of a D meson into the final state of interest is modified from the uncorrelated expectation if the other meson is observed to decay into a CP eigenstate, or indeed any final state that is not flavour specific. These modifications are dependent on the hadronic parameters, and so observables measured in ψ(3770) decays can be used to determine R D and δ S D . This strategy has been pursued on a data set corresponding to an integrated luminosity of 0.81 fb −1 collected by the CLEO-c detector [6][7][8]. Complementary information on the hadronic parameters may be obtained from measurements of D 0 -D 0 oscillations above charm threshold [9]. Here the coherence factor dampens the oscillation signal with respect to that which occurs for the two-body mode D → K − π + . Such a measurement has been performed by the LHCb collaboration for D → K − π + π + π − [10], which, when combined with the results from the CLEO-c data set, yields R K3π = 0.43 +0. 17 −0. 13 [8].
Measurements of CP asymmetries and associated observables in B − → DK − decays using these modes have been performed by LHCb, with data collected during Run 1 of the LHC [11][12][13], and by Belle [14]. The measurements are interpreted using the reported values of the hadronic parameters to obtain constraints on the value of γ and related parameters of the b-hadron decay. In these interpretations the uncertainties associated with the knowledge of the D hadronic parameters are currently smaller than those associated with the finite size of the b-decay samples, but this will no longer be the case with larger Run 2 data set now under analysis at LHCb, and the data sets that both LHCb and Belle II expect to collect over the coming years [15][16][17]. Therefore, more precise measurements of the D hadronic parameters are essential to enable improved knowledge of the angle γ. Such measurements, when obtained at charm threshold, will also be valuable for interpretation of D 0 -D 0 oscillation measurements performed at higher energies.
It is noteworthy that the coherence factor of D → K − π + π + π − is low, which is an attribute that dilutes its sensitivity to γ in a B − → DK − analysis. In order to ameliorate this problem, it is advantageous to partition the phase space of the decay into regions of higher individual coherence, with well-separated values of average strong-phase difference, A binning scheme to achieve this goal has recently been proposed, based on input from D 0 → K − π + π + π − and D 0 → K + π − π − π + amplitude models constructed by LHCb [18]. Interpreting the b-decay properties within these bins requires corresponding binned information for the D hadronic parameters. Initial measurements of these binned parameters has been performed using the CLEO-c data set, but these are rather imprecise [19]. Once more, improved measurements are desirable.
In this paper measurements are reported of the coherence factor and average strongphase differences of D → K − π + π + π − and D → K − π + π 0 decays, integrating over the phase space of both modes. Results are also presented in bins of D → K − π + π + π − phase space. The analysis exploits a data set equivalent to an integrated luminosity of 2.93 fb −1 [20,21] delivered by the BEPCII collider and collected by the BESIII experiment in e + e − collisions at a collision energy of 3773 MeV. The definition of the measured observables and their relationship to the underlying physics parameters is given in Sec. 2. The BESIII detector is described in Sec. 3, together with information concerning the Monte Carlo simulation samples used. The event selection is presented in Sec. 4. The values of the measured observables, the fit to the hadronic parameters, and the assignment of systematic uncertainties are discussed in Sec. 5. Section 6 assesses the impact of the results on measurements of γ in B − → DK − decays, and a summary is given in Sec. 7.

Formalism and measurement strategy
Consider a DD pair produced in e + e − collisions at the ψ(3770) resonance, where one charm meson decays into the signal mode S, either K − π + π + π − or K − π + π 0 , and the other decays into a tagging mode, denoted T . The two mesons are produced in a C-odd eigenstate and their decays are quantum correlated, with a rate given by Note that in general the tagging mode may be a multi-body decay, and therefore has its own coherence factor R T , average amplitude ratio r T D and strong-phase difference δ T D . Throughout the discussion it is assumed that CP violation can be neglected in the charm system. In addition, Eq. 2.1 omits terms of order O(x 2 , y 2 ) ∼ 10 −5 associated with D 0 -D 0 oscillations, where x and y are the usual mixing parameters [22]. These are safe approximations at the current level of experimental precision.
When both S and T are reconstructed the event is said to be double-tagged. Three separate classes of tag are employed in the analysis: CP tags, like-sign tags and K 0 S π + π − tags. These are described for the global measurement, integrated over the full phase space of the signal decay, and then discussed for the case where the D → K − π + π + π − phase space is binned.

CP tags
In the case when the tagging mode is a CP eigenstate, then r T D = 1, R T = 1, and δ T D = 0 or π, and Eq. 2.1 becomes with CP indicating a specific mode with eigenvalue λ (λ = +1 when δ T D = 0 and −1 when δ T D = π). Some self-conjugate decay modes are not pure CP eigenstates, but are known to be predominantly CP even or CP odd. Then Eq. 2.2 continues to apply with the substitution λ = (2F T + − 1), where F T + is the CP -even fraction of the decay. The most striking known example is D → π + π − π 0 , where F πππ o + is in excess of 95% [23,24]. It is useful to define the observables ρ S CP ± , which are the ratios of the number of CPtagged signal events to the number expected in the absence of quantum-correlations: Here N (S|CP ) (N (S|CP )) is the number of decays of channel S (S) tagged with a CP eigenstate after efficiency correction, N DD is the produced number of neutral DD pairs in the sample and B(D 0 → X) is the branching fraction of a D 0 meson into state X. The observable ρ S CP + (ρ S CP − ) corresponds to the case when the CP tag is even (odd).
In relating the squared amplitudes to branching fractions it is necessary to include corrections arising from D 0 -D 0 oscillations. Initially keeping terms to O r 2 D , x 2 , y 2 , the below relations apply: where it is noted that both x and y are of magnitude ∼ r S D /10 [25]. It thus follows that (2.7) This and subsequent expressions are now written to O r 2 D , x, y . In order to combine the information from ρ S CP + and ρ S CP − into a single, CP -invariant, parameter, it is convenient to In practice, the precision of the measurement of ρ S CP ± with many CP eigenstates is limited by the uncertainty on the branching fraction of the tagging mode. This problem can be largely circumvented by re-expressing ρ S CP in a manner which involves N (K ∓ π ± |CP ), the efficiency-corrected yield of D → K ∓ π ± decays tagged by the same CP eigenstate: (2.10) The values of ρ Kπ CP ± can be calculated from Eq. 2.7 with a precision of 0.5%.

Like-sign tags
Consider the case where the tagging mode is also a flavoured final state of the category K − nπ, for example K − π + , K − π + π + π − or K − π + π 0 , and the charge of the kaon in the tagging mode is the same as that of the signal.
In the case where the signal and tagging mode are identical Eq. 2.1 reduces to Again, an observable is constructed that expresses the ratio of the number of efficiencycorrected double-tagged events, N (S|S) and N (S|S), to the expectation without quantum correlations: Since quantum-correlation effects are negligible in opposite-sign double tags, it is experimentally advantageous to determine this observable through the equivalent expression: .
(2.13) Consideration of Eq. 2.1, Eq. 2.11 and the relationship between the squared amplitudes and branching ratios leads to , (2.14) from which it may be seen that this observable has high sensitivity to the coherence factor. When the tag is the two-body mode D → K − π + then R T = 1 and the other parameters r Kπ D and δ Kπ D are well known from measurements of D 0 -D 0 oscillations. The case S = K − π + π + π − and T = K − π + π 0 carries simultaneous information on both multi-body signal modes of interest. The quantum-correlation effects in these double tag can be studied through the observables where T = S. Again, it is beneficial to take advantage of the yields of opposite-sign double tags and to evaluate the equivalent expression It follows from Eq. 2.1, Eq. 2.15 and the branching-ratio relations that 2.3 K 0 S π + π − tags The self-conjugate decay mode D → K 0 S π + π − has been extensively studied at charm threshold, and measurements have been performed by both the CLEO and BESIII Collaborations [26][27][28] of the strong-phase variation over the phase space of the three-body final state.
The Dalitz plot of the decay has axes corresponding to the squared invariant masses The strong-phase difference between symmetric points in the Dalitz plot is given by ∆δ . The bin boundaries are chosen such that each bin spans an equal range in ∆δ K 0 S ππ D (the 'equal-∆δ D binning scheme'), as shown in Fig. 1 where the variation in ∆δ K 0 S ππ D is assumed to follow that predicted by an amplitude model [29]. Measurements performed with quantum-correlated DD pairs determine the actual strong-phase difference within each Dalitz bin. More precisely, what are measured are c i , the cosine, and s i , the sine of the strong-phase difference weighted by the D-decay amplitude in bin i: with an analogous expression for s i . When employing D → K 0 S π + π − as a tag mode it is also necessary to know K i , which is the probability of a single D 0 decay occurring in bin i with no requirement on the decay of the other charm meson in the event: where the sum in the denominator is over all bins. This quantity may be measured in flavour-tagged decays, either at charm threshold or at higher energies (although in the latter case small corrections must be applied to account for the effects of D 0 -D 0 oscillations [30]). Events where one meson decays to K − nπ (K + nπ) and the other to K 0 S π + π − are labelled with a negative (positive) bin number if m 2 − < m 2 + (m 2 − > m 2 + ). With these definitions, it can be shown from Eq. 2.1 that, with uniform acceptance over the Dalitz plot, the yield of double-tagged events is given by where H is a bin-independent normalisation factor. In principle, the channel D → K 0 L π + π − may be employed as a tagging mode in an analogous manner. However, the K i factors for this decay are less well known than those for D → K 0 S π + π − , and limit the sensitivity of this tag. Hence only the D → K 0 S π + π − tag is used in the current analysis.

Binning the D → K
In Ref. [19] it is demonstrated how the sensitivity to γ in B − → DK − , D → K − π + π + π − decays can be enhanced by dividing the phase space of the D decay into bins, each with its own coherence factor and strong-phase difference. A binning scheme is proposed based on the resonance sub-structure of the CF and DCS decays as modelled in the LHCb study reported in Ref. [18]. The phase space is divided into four bins, each of which span a variation in the strong-phase difference between the D 0 andD 0 decay amplitudes. The boundaries are chosen to equalise the product of the integrated CF amplitude and integrated DCS amplitude between bins, which is a metric that is shown to optimise the γ sensitivity. In this binning scheme a small region of phase space is excluded where the invariant mass of π + π − pairs lies close to the mass of the K 0 S meson. This is done in order to exclude background from D → K 0 S K − π + decays, and, according to the amplitude models, rejects 5% of the signal. The phase space of D → K − π + π + π − decays is five-dimensional, and therefore there is no convenient way to visualise the bins. However, the bin in which each decay lies can be unambiguously assigned through knowledge of the four-momenta of the final-state particles.
The observables discussed may be determined within each bin in an identical manner to the global case, giving four measurements each for ρ S CP ± , ρ S F,LS and Y S i . For ρ K3π LS , where the signal and tag decays are both D → K − π + π + π − , ten distinct measurements can be performed.

Detector and simulation samples
The BESIII detector is a magnetic spectrometer [31] located at the Beijing Electron Positron Collider (BEPCII) [32]. The cylindrical core of the BESIII detector consists of a heliumbased multilayer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon-identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% of the 4π solid angle. The charged-particle momentum resolution at 1 GeV/c is 0.5%, and the specific energy loss (dE/dx) resolution is 6% for electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (end-cap) region. The time resolution of the TOF barrel section is 68 ps, while that of the end cap is 110 ps.
Simulated samples produced with the geant4-based [33] Monte Carlo package, which includes the geometric description of the BESIII detector and the detector response, are used to determine the detection efficiencies and to estimate the backgrounds. The simulation includes the beam-energy spread of 0.97 MeV and initial-state radiation (ISR) in the e + e − annihilations modelled with the generator kkmc [34]. The inclusive Monte Carlo samples consist of the production of D 0D0 and D + D − pairs from decays of the ψ(3770), decays of the ψ(3770) to charmonia or light hadrons, the ISR production of the J/ψ and ψ(3686) states, and the continuum processes incorporated in kkmc [34]. The equivalent integrated luminosity of the inclusive Monte Carlo samples is about ten times that of the data. The known decay modes are modelled with evtgen [35,36] using branching fractions taken from the Particle Data Group [37], and the remaining unknown decays from the charmonium states with lundcharm [38,39]. The final-state radiation (FSR) from charged final-state particles is incorporated with the photos package [40]. The signal processes are generated separately taking the spin-matrix elements into account in evtgen. The signal mode D → K − π + π 0 is generated with a resonant sub-structure which matches that of the CF process reported in Ref. [37]; for the decay D → K − π + π + π − the simulated events are re-weighted so that the invariant-mass distributions agree with those produced by the CF model of Ref. [18]. Sample sizes of 200,000 events are simulated for each class of double tag, apart from the case where the tagging mode is D → K 0 S π + π − , for which 400,000 events are produced for each signal channel.

Selection of double-tagged events
In order to determine the observables introduced in Sec. 2, double-tagged samples are collected by reconstructing one charm meson in its decay to a signal mode, and the other charm meson in its decay to a tag mode, which is classed as either flavour, CP , or selfconjugate. Flavour tags can be divided into like sign, where the kaon is of the same sign as the signal mode, or opposite sign. The latter category is used for normalisation purposes. Six categories of CP -even tags are reconstructed, including the quasi-eigenstate mode D → π + π − π 0 , and six CP -odd tags. Events containing both D → K − π + and CP tags are also reconstructed for normalisation purposes. The self-conjugate category contains one tag mode D → K 0 S π + π − . The full list of tags used in the analysis is given in Table 1. Events containing a K 0 L meson are reconstructed using a missing-mass technique, discussed below. All other double-tagged events are fully reconstructed. CP -violation and matterinteraction effects within the neutral-kaon system are not considered because their impact on the measurements is negligible in comparison to the experimental sensitivity. Those mesons that decay within the detector are reconstructed through the modes: K 0 S → π + π − , π 0 → γγ, η → γγ and π + π − π 0 , ω → π + π − π 0 , η → π + π − η and γπ + π − , and φ → K + K − . Table 1. Summary of tag modes selected against the signal decays D → K − π + π + π − and D → K − π + π 0 . The CP tags are also reconstructed against the decay D → K − π + .

Flavour
Like sign Selected charged tracks must satisfy | cos θ| < 0.93, where θ is the polar angle with respect to the beam axis. The distance of closest approach of the track to the interaction point (IP) is required to be less than 10 cm in the beam direction and less than 1 cm in the plane perpendicular to the beam, except for tracks from K 0 S candidates where the closest approach to the IP is required to be within 20 cm along the beam direction. Separation of charged kaons from charged pions is implemented by combining the dE/dx measurement in the MDC and the time-of-flight information from the TOF. This information is used to calculate the probabilities P K and P π for the K and π hypothesis, respectively, and the track is labelled a K (π) candidate if P K > P π (P π > P K ).
Photon candidates are selected from showers deposited in the EMC crystals, with energies larger than 25 MeV in the barrel (| cos θ| < 0.8) and 50 MeV in the end cap (0.86 < | cos θ| < 0.92). In order to suppress fake photons from beam background or electronic noise, the shower clusters are required to be within [0, 700] ns of the start time of the event.
Furthermore, the photon candidates are required to be at least 20 • away from any charged tracks to eliminate fake photons caused by the interactions of hadrons in the EMC.
Candidate K 0 S mesons are reconstructed from pairs of tracks with opposite charge. To improve efficiency no particle-identification requirements are imposed on these tracks. The K 0 S candidate is required to satisfy the flight-significance criterion L/σ L > 2, where L is its flight distance obtained from a fit to the vertex of the track pair, and σ L the uncertainty on this quantity. In addition, a constrained vertex fit is performed for each candidate, and the resulting invariant mass is required to lie within [0.487, 0.511] GeV/c 2 .
In forming π 0 (η) candidates with pairs of photons it is required that the di-photon invariant mass lies within [0.115, 0.150] ([0.480, 0.580]) GeV/c 2 , and that at least one photon candidate is found in the barrel region, where the energy resolution is best. To improve momentum resolution, a kinematic fit is performed, where the reconstructed π 0 (η) mass is constrained to the nominal value [37], and the resulting four-vector is used in the later steps of the analysis. When reconstructing η → π + π − π 0 decays, the invariant mass of the π + π − π 0 combination is required to lie within To suppress combinatorial background, the energy difference, ∆E = E D − √ s/2 is required to be within ±3σ ∆E around the ∆E peak, where σ ∆E is the ∆E resolution and E D is the reconstructed energy of a D candidate in the rest frame of the initial e + e − collision. The resolution varies between decay modes, and the allowed interval in ∆E ranges from To suppress background from cosmic and Bhabha events in the tag modes D → K + K − , π + π − and K − π + , the event is required to have two charged tracks with a TOF time difference less than 5 ns and that neither track is identified as an electron or a muon. As discussed below in Sec. 4.2, a K 0 S veto based on flight distance is applied to D → π + π − π 0 candidates to suppress background from D → K 0 S π 0 , and similarly to D → K − π + π + π − candidates, to reduce contamination from D → K 0 S K − π + decays. If there are multiple double-tagged candidates in one event, the combination with the average reconstructed invariant mass lying closest to the nominal D 0 mass is chosen. Around 10% of selected events contain more than one candidate.
For each fully reconstructed D candidate the beam-constrained mass, M BC , is calculated in order to provide optimal separation between signal and background: where p D is the momentum of the D candidate in the rest frame of the initial e + e − collision. Figures 2 and 3 show M BC distributions of the signal decay for a selection of double tags. The distributions for the other double tags may be found in Appendix A. Double-tagged events where the CP -tag mode involves a K 0 L meson cannot be fully reconstructed. Instead, these events are selected using a missing-mass technique. First the signal mode is reconstructed, and its momentum, p S , is measured in the centre-of-mass frame of the e + e − collision. If more than one candidate is found, that one with the smallest -12 - value of |∆E| is retained. Then the total energy, E X , and momentum, p X , of the charged particles and π 0 candidates not associated with the signal mode are determined. This information allows the missing-mass squared, to be calculated, which is expected to peak at the squared mass of the K 0 L meson for the CP tags under consideration. To suppress contamination from K 0 S → π 0 π 0 decays and other backgrounds, events are rejected that contain surplus π 0 candidates, surplus charged tracks, -13 - any η → γγ candidates, or multiple π 0 candidates that share common showers. Figure 4 shows M 2 miss distributions for those double tags containing a K 0 L meson.

Signal yield determination and consideration of peaking backgrounds
The signal yields for the double tags are determined from an extended unbinned maximumlikelihood fit to the M BC distributions of the tag decay in the case of the fully reconstructed events, and to the M 2 miss distributions for the events containing a K 0 L candidate. (An alternative approach, where a two-dimensional fit is performed to the M BC distributions of both the signal and tag decays in fully reconstructed events, is found to give very similar results, and is further considered when discussing the assignment of systematic uncertainties in Sec. 5.1.2.) The signal shape is modelled from fits to Monte Carlo simulation using the RooKeysPdf one-dimensional kernel estimator [41], and is convolved with a Gaussian function to account for differences between the resolution in data and simulation, and whose width and mean are free parameters in the fit. The combinatorial background of the M BC distribution is described with an ARGUS function [42], and that of the M 2 miss distribution with a RooKeysPdf-determined shape taken from simulation. In general the contribution from background that peaks under the signal region is fixed according to what is found in the Monte Carlo simulation, with a shape determined with RooKeysPdf from the same source. Most of these contributions are at a low level in the fully reconstructed events, summing to typically less than 5%, 10% and 10% in the flavour-tagged, CP -tagged and D → K 0 S π + π − -tagged samples, respectively. However some are significantly larger, or require individual treatment, as is explained below. All fits converge with satisfactory residuals.
The two singly Cabibbo-suppressed decays D → K 0 S K ∓ π ± give rise to the same final states as the signal modes D → K ∓ π ± π ± π ∓ , and lead to significant contamination in the like-sign samples for the global measurement. In the binned measurement these decays are rejected by explicitly excluding the K 0 S region from the phase space that is analysed. In the global analysis a mass veto is undesirable, as it would prevent the measurement being representative of all phase space. Instead an alternative K 0 S veto is applied, in which the two pairs of oppositely charged pions in the decay of the signal candidate are fitted to a common vertex in turn, and their flight significance is determined. Events are rejected when this quantity exceeds a value of two for either combination, after which the net selection efficiency for the background events is O(1%). (This requirement is also imposed in the binned analysis, in addition to the mass veto.) The residual contamination is at a similar level to the signal itself in the affected like-sign samples. Its exact contribution is determined by assuming the measured branching fractions [37] and the efficiencies found in Monte Carlo simulation, and also correcting for quantum-correlation effects not present in the simulation. For the latter calculation the hadronic parameters of the background decay are taken from Ref. [43], and those of the other decays contributing to the double tags from Refs. [8,25]. This background accounts for around 1% of events in the samples containing CP and D → K 0 S π + π − tags, and is negligible for the opposite-sign flavour-tagged events. Another important source of background in the like-sign samples arises from oppositesign events in which both a kaon and pion of opposite charge are misidentified. In order to calibrate the level of this contamination, the rate of misidentification is measured in bins of momentum for a sample of opposite-sign events that is selected with no particleidentification requirements on one of the kaon or pion candidates. This sample includes all the double-tag categories listed in Table 1. The results in data are similar to those found in simulation, as can be seen in Fig. 5. Those differences that are observed are applied as corrections in the simulation. The amount of this background depends on the multiplicity of the final states, because of the momentum dependence of the particle identification. It is approximately double the signal size in the K − π + π 0 vs. K − π + sample, and an order of magnitude smaller than the signal contribution in the K − π + π + π − vs. K − π + π + π − sample.
The decay D → K 0 S π 0 is a dangerous background for D → π + π − π 0 tags, as the two modes have opposite CP eigenvalues. Therefore a K 0 S veto is applied based on the flight significance and identical to that imposed for the D → K − π + π + π − selection. The residual background is corrected for quantum correlations, as is the low level of D → π + π − π 0 contamination in the D → K 0 S π 0 sample. The K 0 S π + π − π 0 final state that is used to construct the CP -odd D → K 0 S ω and D → K 0 S η tags has a non-resonant component that is determined to be at the 20% and 10% level, respectively. Studies reported in Ref. [44] suggest that this background is mildly CP odd. Non-resonant D → K 0 S π + π − π 0 and D → K 0 L π + π − π 0 decays comprise around 30% of the D → K 0 L ω sample, and are expected to be CP neutral. These estimates, obtained from simulation, are corroborated by fits to the π + π − π 0 mass spectrum in data, indicating that any quantum-correlation effects are here negligible. The D → K 0 L π 0 π 0 sample suffers background from D → π 0 π 0 π 0 decays, which generates a structure in the M 2 miss spectrum that peaks at low values. The contribution from this component is a free parameter in the fit, with its shape fixed from simulation studies.
The signal yields for the flavour-tagged and the D → K 0 S π + π − tagged events are given in Table 2. Those for the CP tags may be found in Table 3, together with the selection efficiencies that are required to determine the ρ K3π CP ± and ρ Kππ 0 CP ± observables. Table 2. Signal yields for the flavour-tagged and D → K 0 S π + π − tagged events. The uncertainties are statistical only. Tag In order to improve the resolution of the reconstructed position of each decay in phase space, a kinematic fit is applied to the D → K − π + π + π − candidates with the D 0 mass imposed as a constraint. The events are partitioned into the four bins defined in Ref. [19], and the same procedure is followed as above, to determine the signal yields in each bin. By way of example, the signal yields of D → K − π + π + π − tagged with D → K 0 S π 0 are 497.5 ± 23.2, 415.0 ± 21.0, 400.1 ± 20.2 and 578.6 ± 24.8 for bins 1 to 4, respectively.

Efficiency corrected double-tagged yields in
tagged events are divided into 16 bins of phase space of the tag decay. Again, a mass-constrained fit is applied to improve the resolution of the position of the decay in phase space. Monte Carlo studies indicate that the correct bin is assigned 86% of the time for the positive bins, and 94% of the time for the negative bins. The signal yields are determined in each bin and efficiency corrections are applied. These corrections account for both the relative efficiency variation bin-to-bin, which lies within a range of ±15%, and for migration between bins. The resulting signal yields Y K3π i and Y Kππ 0 i are listed in Table 4. This procedure is repeated for the binned D → K − π + π + π − analysis to measure the signal yields Y K3π,j i presented in Table 5. Tag in bin i of D → K 0 S π + π − phase space. The efficiencies integrated over all bins are 12.9% and 13.6% for D → K − π + π + π − and D → K − π + π 0 decays, respectively.  in bin i of D → K 0 S π + π − phase space (rows) and each bin of D → K − π + π + π − phase space (columns

Measurement of the observables and fit to the hadronic parameters
A global analysis is performed, in which the D → K − π + π + π − and D → K − π + π 0 phase spaces are considered inclusively. In addition, a binned analysis is executed, where the D → K − π + π + π − phase space is partitioned according to the four-bin scheme defined in Ref. [19].
The determination of the CP -tag and opposite-sign observables, and their interpretation in terms of the hadronic parameters, requires knowledge of branching fractions, ratios of branching fractions, and other parameters. The values used, and their sources, are given in Table 6. Table 6. Input parameters used in the determination of the observables and hadronic parameters.

Input parameter
Value Reference  [24] When fitting the signal yields tagged by D → K 0 S π + π − decays in bins of phase space, Y K3π i and Y Kππ 0 i , to the expected distribution of events, expressed in Eq. 2.20, it is necessary to know the strong-phase parameters c i and s i , defined in Eq. 2.18. Values for these parameters are taken from the combined results of measurements performed by the BESIII [27,28] and CLEO collaborations [26]. In addition, knowledge is required of the K i parameters, defined in Eq. 2.19. Here, the most precise source of information comes from models fitted to the flavour-tagged yields of D → K 0 S π + π − decays at the BaBar and Belle experiments. Predictions for K i can be obtained from the models, where is a small correction due to the presence of D 0 -D 0 oscillation effects. Table 7 lists the K i results reported in Ref. [46], from which the K i values used in the current study are derived. These results are considered to be reliable on account of the acceptable fit residuals found in that analysis. The accompanying uncertainties are assigned from the statistical precision of the fit residuals in each bin, added in quadrature to the root-mean-square of the variation of the predictions across this and three other models [29,47,48]. Table 7. The K i parameters as calculated from the model of Ref. [45], equivalent to the fraction of flavour-tagged D → K 0 S π + π − decays falling in each bin of the Dalitz plot, including D 0 -D 0oscillation effects, in the equal-∆δ D binning scheme.

Bin
K The CP -tag observables ρ K3π CP ± and ρ Kππ 0 CP ± are determined for each tag according to Eq. 2.10. They take as input the efficiency corrected CP -tagged signal and D → K − π + yields, listed in Table 3, and the correction factors ρ Kπ CP + = 1.119 ± 0.005 and ρ Kπ CP − = 0.881 ± 0.005, which are calculated from knowledge of the external parameters x, y, r Kπ D and δ Kπ D [25]. The results are given in Tables 8 and 9 and are shown graphically in Fig. 6. The results for each signal mode and class of tag are seen to be compatible and are therefore combined together to give a single, average, result for each observable. Full account is taken of correlations between the systematic uncertainties. Following Eq. 2.8 the parameters ∆ K3π CP and ∆ Kππ 0 CP are also calculated, yielding the values shown in  Table 2. These results are also presented in Table 10. The correlation matrix for all the observables may be found in Appendix B, accounting for both statistical and systematic contributions.
The observables are shown graphically in Fig. 7. The values of ∆ K3π CP and ∆ Kππ 0 CP are incompatible with zero, and several of the like-sign observables are incompatible with unity, indicating the presence of significant effects from quantum correlations. The predictions from the global fit to the BESIII data, discussed below in Sec. 5.1.3, are superimposed.

Assignment of systematic uncertainties
The systematic uncertainties reported in Table 10 are estimated by repeatedly varying each input by a Gaussian distribution with a width set to the uncertainty of that component, and determining the shift in the central value. These shifts are then combined, taking account of any correlations that exist.

Observable
Value Observable Value The like-sign observables are constructed by normalising the like-sign yields by the opposite-sign yields in the same final states, ensuring cancellation for most uncertainties associated with selection efficiencies. However, the resonant sub-structure will in general be different for the like-sign events compared to the opposite-sign sample. The consequence of this difference for D → K − π + π + π − decays is estimated by re-weighting the Monte Carlo sample to the DCS model of LHCb, rather than the CF model [18], and re-evaluating the selection efficiency on K − π + -tagged signal events. A relative difference of 3% is found, which is assigned as an uncertainty to all the like-sign D → K − π + π + π − observables. A similar study performed for D → K − π + π 0 decays, which makes use of the DCS model reported in [50], leads to an uncertainty of 3.5% for the corresponding observables associated with this mode.  Figure 7. Results for the CP -tagged and like-sign global observables. Also shown are the predictions from the fit to these observables and the K 0 S π + π − tags.
Yields are determined from fits to the M BC distribution of the signal decay. In order to assign an uncertainty for this procedure, comparisons are made to the results when a two-dimensional fit is applied to the M BC distributions of both the signal and tag decays, for the example the case where the tag mode is D → K 0 S π 0 . Small differences are observed and uncertainties of 1.5% and 0.8% are assigned to ρ K3π CP ± and ρ Kππ0 CP ± , respectively. No uncertainty is applied for the like-sign observables, as any bias cancels between the like-sign yields and the same-topology opposite-sign yields used for normalisation.
The most important source of background is the decay D → K 0 S K − π + , particularly in the case of the like-sign double tags. The uncertainty on this background contribution is assigned by accounting for the knowledge of its coherence factor and mean strong-phase difference, measured to be 0.70 ± 0.08 and (0.1 ± 15.7) • , respectively [43]. These inputs are required to correct for quantum-correlation in the decay rate. Uncertainties on the branching fraction [45] and the selection efficiency are also considered. The uncertainty on this background component is the dominant systematic uncertainty for the like-sign tags.
A final source of uncertainty on the CP -tag and like-sign observables is the knowledge on the input parameters, which is taken from the measurements reported in Table 6.
When analysing the efficiency-corrected signal yields tagged by D → K 0 S π + π − decays in bins of phase space, Y K3π i and Y Kππ 0 i , it is only necessary to have control of the relative efficiency variation. This variation is less than 15% and taken from Monte Carlo simulation.
It is assumed that any biases on these efficiency corrections are negligible. Small uncertainties are assigned associated with the finite size of the simulation samples. Fits to the expected distribution of events must take account of the uncertainties on the D → K 0 S π + π − strong-phase parameters, found in Ref. [27], and K i parameters, given in Table 7.

Fit to the hadronic parameters
A χ 2 fit is performed to the complete set of CP - (Tables 8 and 9), like-sign (Table 10) and K 0 S π + π − -tagged observables (Table 4) in order to determine the underlying physics parameters: In this fit each individual CP -tag result is entered as a separate measurement. In total, therefore, there are 65 observables and six free fit parameters. The systematic uncertainties from external inputs are not included in the uncertainties of the observables. Instead, the auxiliary parameters δ Kπ D , r Kπ D , x, y, F πππ 0 + , K i , c i and s i are also fitted, with Gaussian constraints introduced into the χ 2 function according to their measured values and covariances. When terms involve the ratio of branching fractions B( to the theoretical predictions. In summary, therefore, the function that is minimised is where χ 2 CP , χ 2 LS and χ 2 K 0 S ππ are the contributions from the CP tags, like-sign tags and K 0 S π + π − tags, respectively, and χ 2 aux contains the constrained contributions from the auxiliary parameters and the ratios of branching fractions. Studies performed with simulated data sets demonstrate that the fit is unbiased and returns correctly estimated uncertainties.
The fit converges with a χ 2 /n.d.f. of 61/59, indicating compatibility between the inputs. Figure 8 shows the Y K3π i and Y Kππ 0 i observables with the fit results superimposed. Also shown are the predictions from the uncorrelated hypothesis, which are directly proportional to the K i values.
The results for the hadronic parameters are given in Table 11 (with the accompanying correlation matrix available in Appendix B), and ∆χ 2 scans in (R K3π , δ K3π D ) and (R Kππ 0 , δ Kππ 0 D ) parameter space are presented in Fig. 9. The measurement of R K3π is in agreement with the value of 0.46 predicted by the LHCb models [18]. These results are also compatible with those obtained from CLEO-c data [8], with an improvement in the 1σ uncertainties for the coherence factors. 3 The region of (R K3π , δ K3π D ) parameter space encompassed by the 2σ and 3σ confidence intervals is significantly more constrained. Combined fits of the BESIII data, together with the CLEO-c and LHCb observables are reported in Appendix C. The fit results for the auxiliary parameters are all compatible with their measured values, and in no case is their precision significantly improved by the fit.
It is instructive to estimate the intrinsic statistical precision of the sample, and the contribution of each source of systematic bias to the overall uncertainty. Therefore, the fit is re-performed many times, with each input varied in turn according to a Gaussian distribution with width set to its assigned uncertainty, and all the other inputs fixed. The width of the distribution of the fitted parameters is taken as an estimate of the uncertainty ) parameter space, showing the ∆χ 2 =2.30, 6.18, 11.83 intervals, which correspond to 68.3%, 95.4% and 99.7% confidence levels in the two-dimensional parameter space. Also shown are the equivalent contours determined from the CLEO-c data [8].
associated with the varying input. The results are presented in Table 12. The most important contributions are seen to come from the finite size of the CP -tagged D → K − π + samples, the uncertainty on the D → K 0 S K − π + background, and on the knowledge of the D → K 0 S π + π − strong-phase and K i parameters. The statistical uncertainty is dominant for all four measurements.

Binned
The binned analysis proceeds in an identical manner to the global case. In this case there are 170 observables and 15 free fit parameters. Because the binning scheme is constructed to exclude D → K 0 S K − π + background there is no uncertainty from this source. In Monte Carlo simulation it is found that around 90% of decays are assigned to the correct bin. A migration matrix, determined from simulation, is used to correct for incorrect assignments. Fits to ensembles of simulated experiments confirm that the results are unbiased and assigned reliable uncertainties. The measured values of the observables are presented in Table 13. The accompanying correlation matrix may be found in Appendix B.
The results for the fit to hadronic parameters are given in Table 11 (with the correlation matrix in Appendix B), and ∆χ 2 scans in (R K3π , δ K3π D ) space are shown in Fig. 10. The fit quality, with χ 2 /n.d.f. = 180/155, is satisfactory. In Appendix C may be found the results for a combined fit to the BESIII and CLEO-c data.
The amplitude models may be used to calculate predictions for the coherence factor in each bin, and the variation in strong-phase between bins [19]. By making use of the measured value of δ K3π D from the global analysis it is then possible to obtain an effective prediction of the average strong-phase difference bin-by-bin, and correlated uncertainty. Following this procedure the predicted values of the coherence factors and strong-phase differences Table 12. Indicative sizes of the systematics uncertainties on the global hadronic parameters. First are listed those components associated with the detector and measurement procedure, and then those associated with the external parameters. Also shown is the statistical uncertainty for comparison.   Table 13. The measured binned observables for D → K − π + π + π − . The first uncertainty is statistical and the second systematic.  6 Impact of the results on the measurement of γ Improved knowledge of the global coherence factors and average strong-phase differences in the processes D → K − π + π + π − and D → K − π + π 0 will be valuable in future studies of B − → DK − decays at LHCb [16,24] and Belle II [17]. However, neither of these channels, when used in isolation, is able to provide a standalone measurement of the Unitarity Triangle angle γ with interesting sensitivity; rather, they bring important constraints when considered as part of the full ensemble of B − → DK − measurements using several D decay modes. These constraints will be tightened as a consequence of the BESIII measurements.

Systematics
Subdividing the phase space of D → K − π + π + π − decays into four bins allows a standalone measurement of γ to be performed [19], and here it is possible to calculate the impact of the measurements reported in this paper. Within each bin of phase space, and considering all possible charge configurations, there are four decay rates that may be measured: The phase difference between these two amplitudes is (δ B − γ), where δ B is a CP -conserving strong phase. These expressions neglect small corrections from D 0 -D 0 oscillations, which can be included in a straightforward manner [51]. It can be seen that the suppressed pair of decays, B ∓ → (K ± π ∓ π ∓ π ± ) D K ∓ , has highest sensitivity to γ, as in this case the interference term involving this parameter enters at leading order. An ensemble of simulated data sets is generated, each containing around 600 suppressed decays. The sample size is roughly equivalent to that which is expected in the Run 1 and Run 2 data sets of LHCb, extrapolating from published results [12]. The distribution of events between B − and B + is simulated according to Eq. 6.1, with the central values of the D hadronic parameters chosen to be the values measured by BESIII, and with the parameters of the B decay and γ = 71 • set to the central values of a recent global average [25]. Each data set is subjected to a χ 2 fit to determine γ, the B-meson decay parameters, and the hadronic parameters of the D decay, where the latter are constrained according to the central values and covariances of the BESIII measurement, as expressed in the likelihood contours. In order to quantify the contributions to the overall fit uncertainties that are induced by the BESIII measurements alone, the exercise is then repeated, but with simulated B-decay data sets that are 100 times larger. Figure 11 shows ∆χ 2 , the change in χ 2 with respect to the minimum that these studies give, plotted as a function of γ. The minimum is centred on the input value. 4 It is seen that an analysis of 600 suppressed B decays, together with the knowledge of the D hadronic parameters reported in this paper, will allow γ to be determined with a precision of +7 −9 • , to which the BESIII measurements contribute an uncertainty of +5 −7 • . Such a result would be only slightly less precise than the current best standalone determination of γ, which comes from an LHCb analysis of B − → DK − , D → K 0 S π + π − and D → K 0 S K + K − decays, and has an uncertainty of ±5 • [52]. Hence it is concluded that the BESIII results on the binned D → K − π + π + π − hadronic parameters are a valuable contribution to improving the knowledge of CP violation in b-hadron decays.  Figure 11. The distribution of ∆χ 2 vs. γ as fitted from two ensembles of simulated data sets, each data set containing 600 or 60,000 suppressed B decays.

Conclusion
Observables have been measured for the decays D → K − π + π + π − and D → K − π + π 0 in a sample of quantum-correlated DD pairs. Many of these observables exhibit significant deviations from the values expected in the absence of quantum correlations, and the ensemble of measurements allow for the coherence factor and average strong-phase difference of each decay to be determined. In the case when the analysis integrates over the phase space of both decays, the results are • , The measurements of the coherence factor are more precise than the existing world-average values, and the allowed region in the (R K3π , δ K3π D ) and (R Kππ 0 , δ Kππ 0 D ) planes is more restricted than is the case for the observables measured with CLEO-c data [8]. This improved knowledge will be valuable when these channels are included with other D-meson decay modes in studies of B − → DK − decays at LHCb and Belle II, and will enable a more precise determination of the angle γ of the Unitarity Triangle. It will also be useful in interpreting D 0 -D 0 oscillation measurements performed with these multibody decays.
The analysis has been re-performed in four bins of phase space of the D → K − π + π + π − decay, to yield bin-by-bin values of the coherence factor and average strong-phase difference. A study of B − → DK − decays that makes use of such a sub-division is inherently more sensitive, and can provide a standalone measurement of γ that is expected to be among the most precise available using a single D-decay mode. The BESIII results will contribute an uncertainty of around 6 • to this measurement, which is less than that estimated to arise from the size of currently available B-meson decay samples. Updated analyses with the larger data sets that BESIII expects to collect at the ψ(3770) resonance in the coming years will allow this uncertainty to be significantly decreased [53].   Figure 12. M BC distributions for decays tagged by the modes D → K 0 S π + π − , D → K 0 S π 0 and D → K 0 S π 0 π 0 . The points with error bars are data; the red line indicates the total fit; the dashed purple line shows the peaking-background contributions; the combinatorial-background contribution is at too low a level to be visible.  Figure 13. M BC distributions for decays tagged by the modes D → K 0 S ω(π + π − π 0 ), D → K 0 S η(γγ) and D → K 0 S η(π + π − π 0 ). The points with error bars are data; the red line indicates the total fit; the dashed purple line shows the peaking-background contributions; the combinatorial-background contribution is at too low a level to be visible.
The points with error bars are data; the dashed purple line shows the peaking-background contributions; the combinatorial-background contribution is at too low a level to be visible. Table 14 presents the correlation matrix for the observables in the global analysis, and  Table 15 gives the correlation matrix for the hadronic parameters determined in this study. Tables 16 and 17 show the corresponding matrices for the binned analysis.

C Combination with CLEO-c and LHCb results
A fit is performed to the global observables determined by BESIII and those from the CLEO-c data, reported in Refs. [7,8]. It is assumed that the correlations between the two sets of measurements are negligible. The χ 2 /n.d.f. of the fit is 100/98. The results are presented in Table 18, with the accompanying correlation matrix given in Table 19, and ∆χ 2 scans are shown in Fig. 15. A second fit, with a χ 2 /n.d.f. of 104/101, includes the constraints from the LHCb study of D 0 -D 0 oscillations [10]. The results are shown in Table 18 and Fig. 16, with the correlation matrix in Table 20.   Table 21 and Fig. 17. The correlation matrix is given in Table 22.