Quantum-correlated measurements of $D\to K^{0}_{\rm S}\pi^{+}\pi^{-}\pi^{0}$ decays and consequences for the determination of the CKM angle $\gamma$

We perform quantum-correlated measurements of the decay $D \to K^{0}_{\rm S}\pi^{+}\pi^{-}\pi^{0}$ using a data sample corresponding to an integrated luminosity of 0.82 fb$^{-1}$ collected at the $\psi(3770)$ resonance by the CLEO-c detector. The value of the $CP$-even fraction $F_{+}$ is determined to be 0.238 $\pm$ 0.012 $\pm$ 0.012. The strong-phase differences are also measured in different regions of $K^{0}_{\rm S}\pi^{+}\pi^{-}\pi^{0}$ phase space by binning around the intermediate resonances present. The potential sensitivity of the results for determining the CKM angle $\gamma$ from $B^{\pm} \to D(K^{0}_{\rm S}\pi^{+}\pi^{-}\pi^{0})K^{\pm}$ decays is also discussed.


Introduction
The CKM [1,2] angle γ ≡ arg(−V ud V * ub /V cd V * cb ), sometimes denoted as φ 3 , can be measured using decays B ± → DK ± with D being a neutral charm meson reconstructed in a final state common to both D 0 andD 0 decays. The current uncertainty on γ is significantly larger than that of the Standard Model (SM) prediction, which is calculated from the assumption of unitarity and measurements of other parameters in the CKM matrix [3]. This is due to the small branching fraction of decays sensitive to γ. A more precise measurement of γ is crucial for testing the SM description of CP violation and probing for new physics effects. The data collected at detectors such as BABAR, Belle, LHCb or the future Belle II experiment can be used to determine γ. The D decay final states so far studied include those that are either CP -eigenstates, flavour specific or self-conjugate. The statistical uncertainty on γ can be reduced if information from additional D meson final states is included, which in practice means new three and four-body decay modes. However, the use of multibody final states requires knowledge of the strong-phase difference between the D 0 andD 0 that varies over the phase space. Determining the fractional CP content of the D meson final state is also helpful. The required information can be obtained by studying quantum-correlated DD mesons produced in e + e − collisions at an energy corresponding to the mass of the ψ(3770).
Here, we present the first quantum-correlated analysis of the decay D → K 0 S π + π − π 0 , which has a branching fraction of 5.2% [3], which is large compared to that of other multibody final states. The study is made with data collected by the CLEO-c detector, corresponding to an integrated luminosity of 0.82 fb −1 . We determine the CP -even fraction F + of the decay which makes it potentially useful in a quasi-GLW [4,5] analysis along with other CP eigenstates [6]. Furthermore, this multibody self-conjugate decay occurs via many intermediate resonances, such as K 0 S ω and K * ± ρ ∓ , hence if the strong-phase difference variation over the phase space is known, a GGSZ-style [7,8] analysis to determine γ from this final state alone is possible. Hence measurements of the relative strong-phase are performed in localised regions of phase space. A similar study has already been performed for a four-body D decay to π + π − π + π − [9].
The remainder of the paper is arranged as follows. Section 2 describes the quantumcorrelated D mesons produced at CLEO-c. The relations used to determine F + as well as the strong-phase differences are also discussed here. The data set and event selection criteria are explained in section 3. The F + and strong-phase difference calculations and results are presented in section 4 and section 5, respectively. The impact of these results on the determination of γ, illustrated with simulated experiments using the expected yield from the data set collected by the Belle II experiment, is discussed in section 6. Section 7 gives the conclusions.

Quantum-correlated D mesons
Decays of the vector meson ψ(3770) produce pairs of D mesons in a P -wave state and hence the wave function for the decay is antisymmetric. It is possible to show that the decay rate is maximum when both the D mesons decay to states of opposite CP eigenvalue and zero when both have same CP eigenvalue [6].

CP -even fraction F +
The CP -even fraction F f + of a D decay to final-state f can be determined from samples in which both D mesons are fully reconstructed; these we term "double-tagged" events. In the case that one decay is reconstructed in final state f and the other to a CP eigenstate g, the double-tagged yield in the integrated phase space can be written in terms of F + , the CP eigenvalue λ g CP and the branching fractions for the states D → f and D → g, B(f ) and B(g) as where N is the number of D 0D0 pairs and (f |g) is the reconstruction efficiency. The single-tagged yield, where only one of the D meson decays is reconstructed, is given by Then F + can be defined as for CP -odd and CP -even g modes, respectively. We can also use some multibody modes g with already known CP -even fraction F g + , to determine F f + [10]: where N g is the ratio of double-tagged to single-tagged yields. The tagging mode g can also be studied in bins of its own phase space, e.g. in the case of K 0 S π + π − or K 0 L π + π − , such that the yield in the i th bin of the K 0 S,L π + π − Dalitz plot is given by are the fractional rate of the D 0 (D 0 ) decays to K 0 S,L π + π − and the cosine of the average strong-phase difference, respectively, in the i th bin [10]. Here, h K 0 S,L π + π − is a normalization factor. Events where both the D mesons decay to the same final state f also provide useful information about F + . The double-tagged yield in that case is given by [10] (2.6)

Strong-phase difference
The sine and cosine of the amplitude weighted averages of the strong-phase difference between D 0 andD 0 , are represented as c i and s i , in the decay D → K 0 S π + π − π 0 , similar to that in D → K 0 S π + π − given in section 2.1. The amplitude for D 0 → f can be written as A(D 0 → f (x)) ≡ a x e iθx , where x is some point in the decay phase space and θ x is the strongphase of the decay that conserves CP symmetry. Similarly forD 0 → f , the amplitude is Here,x is a point in phase space obtained by applying a CP transformation to the final state system at x. We can define the strong-phase difference between D 0 andD 0 as ∆θ x = θ x − θx. So, c i and s i are defined as and where T i = x∈i a 2 x dx andT i = x∈i a 2 x dx. Because of the quantum correlation, the decays of D → K 0 S π + π − π 0 recoiling against CP and quasi-CP eigenstates and other self-conjugate states as tag modes provide direct sensitivity to c i and s i . The double-tagged yields for a CP tag can be written as Type Modes CP -even where h CP is a normalization constant and i represents a particular region of the decay phase space of f . For a quasi-CP tag, the c i sensitive term is scaled by (2F + − 1) rather than 1. For the tag modes K 0 S,L π + π − [11,12], the double-tagged yield is (2.10) Here j is a particular region of the decay phase space of K 0 S,L π + π − . If both the D meson final states are the same, then where h f is a normalization constant.

Data set and event selection
A data sample consisting of DD pairs coming from the ψ(3770) resonance collected by the CLEO-c detector at the CESR e + e − collider is used in this analysis. This corresponds to an integrated luminosity of 0.82 fb −1 . A detailed description of the CLEO-c detector is given in refs. [13][14][15][16]. Monte Carlo (MC) simulations of signal events are used to estimate selection efficiencies. Generic samples of DD MC events having twenty times the integrated luminosity of the data set are used to determine the background contributions. The EvtGen [17] package is used to generate the decays and the detector response is modelled with Geant [18]. The final-state radiation effects associated with charged particles are simulated with PHOTOS [19]. One of the D mesons is reconstructed in the final state of interest, K 0 S π + π − π 0 , and the other to one of the different tag states given in table 1. All tracks and showers associated with both the D mesons are reconstructed; the selection criteria for the tag modes are identical to those presented in ref. [6].
Hadronic modes not involving a K 0 L or ν are fully reconstructed using the kinematic variables beam-constrained mass (m bc ) and the beam-energy difference (∆E), which are where E beam is the beam energy and P D and E D are the summed momenta and energy of the D daughter particles, respectively. A kinematic fit is performed to constrain the final state particles to the D meson invariant mass. This fit improves the momentum resolution of the D daughter particles. For a correctly reconstructed D meson, m bc and ∆E peak at the nominal D mass [3] and zero, respectively. No peaking structure is observed for combinatorial backgrounds. The double-tagged yield is calculated by counting the events in the signal and sideband regions of m bc with different selection criteria on ∆E for various modes as in ref. [6]. For the mode K 0 S π + π − π 0 , not previously analysed, we impose −0.025 GeV < ∆E < 0.025 GeV. Figure 1 shows the ∆E distribution for single-tagged K 0 S π + π − π 0 candidates. Double-tagged events containing two K 0 S π + π − π 0 decays are also reconstructed in a similar fashion. The m bc distributions for D → K 0 S π + π − π 0 decays tagged with CP eigenstates, not involving a K 0 L , and K 0 S π + π − are shown in figure 2. An example of the two-dimensional m bc plane distribution for K 0 S π + π − π 0 vs K 0 S π 0 candidates is shown in figure 3.
Modes with a K 0 L meson in the final state cannot be reconstructed fully because the K 0 L escapes the detector before leaving any useful signature. Hence a missing-mass squared technique [20] is used to reconstruct those events. The m 2 miss is calculated as where E miss is the missing energy and P miss is the magnitude of the missing three-momentum in the event. For a correctly reconstructed event, m 2 miss peaks near the square of the K 0 L mass [3]. The double-tagged yields are estimated from the signal and sideband regions of   Figure 3. Two dimensional m bc plane distribution for K 0 S π + π − π 0 tagged with K 0 S π 0 decays. The red square box indicates the signal region and the remaining boxes show the various sideband regions that are used to determine the combinatorial background contribution. m 2 miss distribution. Similarly, semileptonic decays involving a neutrino are reconstructed by considering the quantity which peaks near zero for a correctly selected event. The yield in this category is estimated by looking at the signal and sideband regions of U miss distribution. The m 2 miss distributions for D → K 0 S π + π − π 0 decays tagged by CP -even states involving a K 0 L meson and K 0 L π + π − along with the U miss distribution for K ± e ∓ ν e tag are shown in figure 4 and 5, respectively. The tag-side Dalitz plot distributions for K 0 S,L π + π − are shown in figure 6. In events that contain more than one reconstructed pair of D meson decays, the candidate with average m bc of both the D mesons closest to the nominal mass of D is chosen [3]. A K 0 S veto is applied for the final state π + π − π 0 to eliminate K 0 S (π + π − )π 0 background tags. As MC does not simulate quantum correlations, a correction is applied to obtain the correct amount of contamination from this source. The peaking background estimated from MC samples for modes not involving a K 0 L meson constitute a maximum of 2.5% of the  selected events. For K 0 L π 0 , K 0 L ω and K 0 L π + π − final states, there is contamination from the corresponding K 0 S modes in the signal region. Again, the raw yields found in the MC are adjusted to account for quantum correlations. The mode K ± e ∓ ν contains no significant peaking background contributions. The background subtracted double-tagged yields for each of the modes are given in table 2.
The single-tagged yields for the CP and quasi-CP modes are taken from ref. [6] as the selection criteria applied are the same. The single-tagged yield for K 0 S π + π − π 0 is obtained from a fit to the m bc distribution and found to be 54,949 ± 781. The signal component is modelled with an asymmetric Gaussian and a sum of two Gaussian probability density functions (PDF) with common mean and the background component is fitted with Argus [21], Crystal Ball [22] and Gaussian PDFs. The latter two PDFs in the background fit are for the small peaking component arising from π + π − π + π − π 0 and K 0 S K 0 S π 0 . The m bc The axis labels m 2 ± represents the invariant mass squares of K 0 S,L π ± pairs. Type Mode Yield Quasi-CP π + π − π 0 428.8 ± 21.7 Table 2. Background subtracted signal yields of D → K 0 S π + π − π 0 for different tag modes.
distribution and fit are shown in figure 7.

Measurement of F +
We perform F + measurements with different tag modes using the relations presented in section 2.1. The results are discussed in the following subsections.  Figure 7. m bc distribution for single-tagged D → K 0 S π + π − π 0 decays. The black points are data, the solid blue curve is the total fit and the dashed red and blue curves are signal and background fit components, respectively.

CP and quasi-CP tags method
The double-tagged yields involving a CP eigenstate tag are used to obtain N + and N − . The dependence on branching fraction and reconstruction efficiency is removed by the normalization with single-tagged yields. The possible effect of DD mixing is eliminated by applying the correction factor for single-tagged yields as S = S meas /(1 − λ CP y D ), where y D = (0.69 ± 0.06)% is the D-mixing parameter [23]. The N + and N − values are shown in figure 8. It can be seen that there is consistency among the values obtained from different modes. From these results, a value of F + = 0.240 ± 0.018 ± 0.011 is calculated using eq. (2.3). The uncertainties are statistical and systematic, respectively. This value indicates that the mode K 0 S π + π − π 0 is significantly CP -odd. The dominant systematic uncertainty comes from the determination of the single-tagged yields: in particular the fit shapes, branching fraction and reconstruction efficiency values used for K 0 L modes. These uncertainties are estimated in an identical manner to that described in ref. [6]. Now, using the quasi-CP tag π + π − π 0 , whose F + value is 0.973 ± 0.017 [10], the CPeven fraction for K 0 S π + π − π 0 is calculated using eq. (2.4). The result obtained with this quasi-CP mode is 0.244 ± 0.020 ± 0.007.
The double-tagged decays with K 0 S π + π − and K 0 L π + π − are analysed by dividing the Dalitz plot of the tag mode into eight pairs of symmetric bins as in ref. [24] according to the amplitude model described in ref. [25]. The symmetric bins are folded across the line m 2 + = m 2 − to make a total of eight bins. The double-tagged yield in each of the folded bins is related to F + as given in eq. (2.5). Therefore F + can be extracted from a combined log-likelihood fit to the yields. The reconstruction efficiency in each bin is obtained from simulated signal samples and a correction is applied to the yields to account for the variation of efficiency across the bins, which varies by typically 3%, bin-to-bin. Table 3 shows the background subtracted efficiency corrected yields in each of the eight bins.
A log-likelihood fit is performed with the input yields following the form of eq. (2.5) with the CP -even fraction and overall normalization as fit parameters. The uncertainty on the K i ,K i , c i and c i input parameters are added as Gaussian constraints in the fit. The fit is performed separately for K 0 S π + π − and K 0 L π + π − and then for both the tags combined. All the fits have good quality and the results are presented in table 4. The measured and predicted yields in each bin are given in figure 9 for both tags.
K 0 S,L π + π − 0.255 ± 0.029 1.42 Table 4. F + results for the mode K 0 S π + π − π 0 from the tags K 0 S π + π − and K 0 L π + π − . The row K 0 S,L π + π − indicates that the combined fit includes both the samples. The fit quality metric χ 2 /DoF is also shown, where DoF stands for the number of degrees of freedom.  Figure 9. Predicted and measured yields for (a) K 0 S π + π − and (b) K 0 L π + π − in each bin obtained from the combined fit of both the modes. The histogram shows the predicted values from the fit, points show the measured values, the dashed line corresponds to F + = 0 and the dotted line shows F + = 1.
There is a two standard deviation difference between the results from each of the tags alone, however the combined result agrees with F + from the other tag methods. The nonuniform acceptance of the K 0 S,L π + π − Dalitz plane is studied by varying the efficiency by 3%. The resulting change of +0.007 −0.008 in F + is assigned as the systematic uncertainty related to this source.

K 0
S π + π − π 0 self-tags method and combined result The double-tagged decays in which both the D mesons decay to the same final state of K 0 S π + π − π 0 can also give information about F + following the relation given in eq. (2.6). The value is obtained to be 0.226 ± 0.019 ± 0.004. Here, the systematic uncertainty arises from the uncertainty on external input values used in the calculation such as the number of DD pairs and the branching fraction of the decay.
The value of F + from all these above methods are given in table 5. They are consistent with each other and the combined result obtained via weighted averaging is 0.238 ± 0.012 ± 0.003, where the correlation due to the use of N + for CP tags as well as the π + π − π 0 tag is taken into account.

Method
F + CP tags 0.240 ± 0.018 ± 0.011 quasi-CP tag 0.244 ± 0.020 ± 0.007 We need to consider another source of systematic uncertainty common to all methods: the non-uniform acceptance across the phase space of D → K 0 S π + π − π 0 , which will bias the result with respect to the flat acceptance case. We estimate the acceptance systematic uncertainty by calculating F + from the c i strong-phase difference results given in section 5, which have bin-wise efficiency corrections. The value of F + is related to c i by The same data are used, so any difference can be attributed to the absence of acceptance corrections in the inclusive method. The obtained result is 0.226 ± 0.020. There is a one standard deviation difference between the value obtained from eq. (4.1) and the averaged unbinned F + result. The difference, 0.012, is taken as the systematic uncertainty from this source. Including this uncertainty the combined result becomes 0.238 ± 0.012 ± 0.012.

Determination of c i and s i
The c i and s i values are extracted by looking at the same D → K 0 S π + π − π 0 data in bins of phase space. The decay phase space is five-dimensional, hence there is no trivial symmetry to define the bins as in the case for three-body decays. Furthermore, a proper optimization is impossible due to the lack of an amplitude model. Therefore, the bins are constructed around the most significant intermediate resonances present in the decay. A nine-bin scheme is defined around the intermediate resonances such as the ω, K * and ρ. The kinematic regions of the bins are given in table 6 along with the fractions of flavour-tagged D 0 and D 0 decays in each of them, which are determined from the semileptonic D → K ± e ∓ ν e double-tagged events; the relevant kinematic distributions are shown in figure 10. The bins are exclusive and the cuts are applied sequentially in the order of the bin number. We also note that increasing the number of bins, which would result in better sensitivity to γ, led to instabilities in the fit due to the large number of null bins. MC studies led to robust results for the nine bins. We do not use K ± π ∓ , K ± π ∓ π 0 or K ± π ∓ π ± π ∓ as flavour tags because the corrections from the Cabibbo-favoured and doubly-Cabibbo-suppressed amplitudes cannot be calculated in the absence of an amplitude model for  Table 6. Specifications of the nine exclusive bins of D → K 0 S π + π − π 0 phase space along with the fraction of flavour-tagged D 0 andD 0 events in each of them. m L and m U are the lower and upper limits, respectively, of the invariant masses in each region.
Due to the finite resolution of the detector, reconstructed decays may migrate to other bins in phase space. This effect is studied by looking at simulated signal events and a 9 × 9 migration matrix M , is calculated. Each of its elements gives the ratio of the number of events reconstructed to those generated in a bin. There is a significant loss of 20% from bin 1 due to the ω resonance having a narrow decay width. The double-tagged yields (Y ) for each mode are corrected for the migration effects as Y i = Σ j M ij Y j , where i and j run from 1 to 9. The K i andK i values given in table 6 are also obtained after the correction applied as K i = Σ j M −1 ij K j . The background subtracted and migration corrected double-tagged yields for CP tags, quasi-CP tag and other self-conjugate modes are obtained in each of the bins. These inputs are used in a Poissonian log-likelihood fit with c i and s i values as fit parameters. The fit assumes that the data follow eqs. (2.9) -(2.11). The CP and quasi-CP tags provide sensitivity only to c i values. The tags K 0 S π + π − and K 0 L π + π − give sensitivity to both c i and s i values. The already measured strong-phase parameters for D → K 0 S,L π + π − are used as inputs in the fit. As the binning scheme for the signal mode is not symmetric, it is no longer possible to exploit the symmetry of the tagging decay. Therefore, there are sixteen bins in the tag-side and nine bins in the signal-side. The sample of doubly-tagged K 0 S π + π − π 0 events is also useful in providing information on s i values. For such events, there are nine bins each in the signal and tag side.
is chosen for one CP tag, K + K − , and all the other normalizations for events not involving K 0 L modes are defined as S(tag) S(K + K − ) h CP in the fit, where S represents the single-tagged yield. The nature of the symmetry within the bins leads to certain constraints that can be imposed in the fit. Bins 1, 6 and 9 are CP self-conjugate, which implies s i = 0. The bins 2 and 3, 4 and 5, and 7 and 8 are each CP -conjugate pairs, which imposes relations between their s i values. We have: s 1 = 0, s 6 = 0, s 9 = 0; (5.1) 3) In the fit, we constrain s 3 , s 5 and s 8 using eqs. (5.2)-(5.4) along with fixing s 1 , s 6 and s 9 to zero.

Systematic uncertainties
Several The limited statistics of the MC sample used to determine the migration matrix can cause variations in phase-space acceptance and biases in the results. The elements of migration matrix are smeared by +1% and −1% independently to account for this possible bias. The resulting change in c i and s i are assigned as systematic uncertainty. This 1% deviation is large enough to take care of the effect of choosing a wrong candidate during multiple candidate selection. The single-tagged yields used in the normalization of the fit are fluctuated independently to +1σ and −1σ, where σ is the statistical uncertainty on the yield and the change in c i and s i values are taken as systematic uncertainty.
The uncertainty on K i andK i values is taken as a Gaussian constraint in the fit and hence no need to assign a systematic uncertainty for this. We investigate the change in migration matrix due to momentum resolution. Bin 1, which has the largest migration, hence the largest sensitivity to any data-MC discrepancy is chosen for the study. The π + π − π 0 invariant mass resolution in data and MC are found to be 5.319 ± 0.064 MeV/c 2 and 4.928 ± 0.003 MeV/c 2 , respectively. The invariant mass in data is smeared by the quadrature difference of these two resolutions 2.019 MeV/c 2 and the migration matrix is recalculated. The effect is minimal and hence we do not assign a systematic uncertainty for this in bin 1 or in any other bins.
The multiplicity distribution for each tag shows good agreement between data and MC. The efficiency of the selection of the best candidate in an event is ≥ 83% in each case. So the metric choosing the best candidate in an event does not introduce a bias. A summary of the systematic uncertainty evaluation is given in table 7 and 8. The systematic uncertainties are small compared with the statistical errors.
The final results of the c i and s i values are given in table 9 and displayed graphically in figure 11. The statistical and systematic correlation coefficients between c i and s i values are given in table 10 and table 11   6 Estimation of γ sensitivity with B ± → D(K 0 S π + π − π 0 )K ± In order to estimate the impact of these results on a future γ measurement using B ± → D(K 0 S π + π − π 0 )K ± decays, we perform a simulation study based on the expected yield of this mode in the Belle data sample (≈ 1 ab −1 ). The Belle sample of B ± → D(K 0 S π + π − )K ± [26] has ≈ 1200 events. Assuming that increase in branching fraction for K 0 S π + π − π 0 compared to K 0 S π + π − is compensated by the loss of efficiency due to a π 0 in the final state [27,28], we expect a similar yield for B ± → D(K 0 S π + π − π 0 )K ± . Then 60,000 is the yield extrapolated to the 50 ab −1 data set anticipated at Belle II. The γ sensitivity is estimated in a GGSZ [7,8] framework. We run 1000 pseudo experiments with c i , s i , K i , andK i values as inputs with each experiment consisting of ≈ 60,000 events. The input values of γ and the hadronic parameters r B and δ B are taken from ref. [29]. The estimated uncertainty on γ is σ γ = 4.4 • (see figure 12). This sensitivity is very promising and only a factor two worse than that 0.00 Table 9.  Figure 11. c i and s i values in each bin. The black and red error bars represent statistical and systematic uncertainties, respectively. anticipated from studying B ± → D(K 0 S π + π − )K ± [30] decays.

Conclusions
Improving the knowledge of the CKM angle γ is an important goal in flavour physics. This can be achieved by harnessing new D decay modes for the measurements of CP asymmetries in B ± → DK ± . We present the first measurement of the CP -even fraction F + for the decay D → K 0 S π + π − π 0 which gives F + = 0.238 ± 0.012 ± 0.012. The F + measurement can be used in a quasi-GLW analysis in which there is no binning of the D → K S π + π − π 0 phase space, although this does not provide single-mode sensitivity to γ. In addition, the    Table 11. Systematic correlation coefficients between c i and s i values.
measurements of amplitude weighted averages of the cosine and sine of the strong-phase difference between D 0 andD 0 decaying to the self-conjugate final state of K 0 S π + π − π 0 have been performed. This is done in nine regions of the decay phase space binned according to the intermediate resonances present. These results allow a model-independent GGSZ estimation of γ from this mode alone. It is estimated that a single-mode uncertainty on γ of σ γ = 4.4 • is achievable with a 50 ab −1 sample of data at Belle II experiment. This could be improved with optimized c i and s i values provided a proper amplitude model is available and a finer binning using a larger sample of quantum-correlated data from the BESIII experiment.