Studies of the resonance structure in $D^0\to K^0_S K^{\pm}\pi^{\mp}$ decays

Amplitude models are applied to studies of resonance structure in $D^0\to K^0_S K^- \pi^+$ and $D^0\to K^0_S K^+ \pi^-$ decays using $pp$ collision data corresponding to an integrated luminosity of $3.0\,\mathrm{fb}^{-1}$ collected by the LHCb experiment. Relative magnitude and phase information is determined, and coherence factors and related observables are computed for both the whole phase space and a restricted region of $100\,\mathrm{MeV/}c^2$ around the $K^{*}(892)^{\pm}$ resonance. Two formulations for the $K\pi$ $S$-wave are used, both of which give a good description of the data. The ratio of branching fractions $\mathcal{B}(D^0\to K^0_S K^+ \pi^-)/\mathcal{B}(D^0\to K^0_S K^- \pi^+)$ is measured to be $0.655\pm0.004\,(\textrm{stat})\pm0.006\,(\textrm{syst})$ over the full phase space and $0.370\pm0.003\, (\textrm{stat})\pm0.012\,(\textrm{syst})$ in the restricted region. A search for $CP$ violation is performed using the amplitude models and no significant effect is found. Predictions from $SU(3)$ flavor symmetry for $K^{*}(892)K$ amplitudes of different charges are compared with the amplitude model results.


Introduction
A large variety of physics can be accessed by studying the decays 1 D 0 → K 0 S K − π + and D 0 → K 0 S K + π − . Analysis of the relative amplitudes of intermediate resonances contributing to these decays can help in understanding the behavior of the strong interaction at low energies. These modes are also of interest for improving knowledge of the Cabibbo-Kobayashi-Maskawa (CKM) [1,2] matrix, and CP -violation measurements and mixing studies in the D 0 -D 0 system. Both modes are singly Cabibbo-suppressed (SCS), with the K 0 S K − π + final state favored by approximately ×1.7 with respect to its K 0 S K + π − counterpart [3]. The main classes of Feynman diagrams, and the sub-decays to which they contribute, are shown in Fig. 1.
Flavor symmetries are an important phenomenological tool in the study of hadronic decays, and the presence of both charged and neutral K * resonances in each D 0 → K 0 S K ± π ∓ mode allows several tests of SU(3) flavor symmetry to be carried out [4,5]. The K 0 S K ± π ∓ final states also provide opportunities to study the incompletely understood Kπ S-wave systems [6], and to probe several resonances in the K 0 S K ± decay channels that are poorly established.
An important goal of flavor physics is to make a precise determination of the CKM unitarity-triangle angle γ ≡ arg(−V ud V * ub /V cd V * cb ). Information on this parameter 2 can be obtained by studying CP -violating observables in the decays B − → ( ) D 0 K − , where the D 0 and D 0 are reconstructed in a set of common final states [7,8], such as the modes D 0 → K 0 S K − π + and D 0 → K 0 S K + π − [9]. Optimum statistical power is achieved by studying the dependence of the CP asymmetry on where in three-body phase space the D-meson decay occurs, provided that the decay amplitude from the intermediate resonances is sufficiently well described. Alternatively, an inclusive analysis may be pursued, as in Ref. [10], with a 'coherence factor' [11] parameterizing the net effect of these resonances. The coherence factor of these decays has been measured by the CLEO collaboration using quantum-correlated D 0 decays at the open-charm threshold [12], but it may also be calculated from knowledge of the contributing resonances. In both cases, therefore, it is valuable to be able to model the variation of the magnitude and phase of the D 0 -decay amplitudes across phase space.
The search for CP violation in the charm system is motivated by the fact that several theories of physics beyond the Standard Model (SM) predict enhancements above the very small effects expected in the SM [13][14][15]. Singly Cabibbo-suppressed decays provide a promising laboratory in which to perform this search for direct CP violation because of the significant role that loop diagrams play in these processes [16]. Multi-body SCS decays, such as D 0 → K 0 S K − π + and D 0 → K 0 S K + π − , have in addition the attractive feature that the interfering resonances may lead to CP violation in local regions of phase space, again motivating a good understanding of the resonant substructure. The same modes may also be exploited to perform a D 0 -D 0 mixing measurement, or to probe indirect CP violation, either through a time-dependent measurement of the evolution of the phase 1 The inclusion of charge-conjugate processes is implied, except in the definition of CP asymmetries. 2 Another notation, φ 3 ≡ γ, exists in the literature.
space of the decays, or the inclusive K 0 S K − π + and K 0 S K + π − final states [17]. In this paper time-integrated amplitude models of these decays are constructed and used to test SU(3) flavor symmetry predictions, search for local CP violation, and compute coherence factors and associated parameters. In addition, a precise measurement is performed of the ratio of branching fractions of the two decays. The data sample is obtained from pp collisions corresponding to an integrated luminosity of 3.0 fb −1 collected by the LHCb detector [18,19] during 2011 and 2012 at center-of-mass energies √ s = 7 TeV and 8 TeV, respectively. The sample contains around one hundred times more signal decays than were analyzed in a previous amplitude study of the same modes performed by the CLEO collaboration [12].
The paper is organized as follows. In Sect. 2, the detector, data and simulation samples are described, and in Sect. 3 the signal selection and backgrounds are discussed. The analysis formalism, including the definition of the coherence factor, is presented in Sect. 4. The method for choosing the composition of the amplitude models, fit results and their systematic uncertainties are described in Sect. 5. The ratio of branching fractions, coherence factors, SU(3) flavor symmetry tests and CP violation search results are presented in Sect. 6. Finally, conclusions are drawn in Sect. 7.

Detector and simulation
The LHCb detector is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector surrounding the pp interaction region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. The tracking system provides a measurement of momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary pp interaction vertex (PV), the impact parameter, is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov (RICH) detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter.
The trigger [20] consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, in which all charged particles with p T > 500 (300) MeV/c are reconstructed for 2011 (2012) data. At the hardware trigger stage, events are required to have a muon with high p T or a hadron, photon or electron with high transverse energy in the calorimeters. For hadrons, the transverse energy threshold is 3.5 GeV. Two software trigger selections are combined for this analysis. The first reconstructs the decay chain D * (2010) + → D 0 π + slow with D 0 → h + h − X, where h ± represents a pion or a kaon and X refers to any number of additional particles. The charged pion originating in the D * (2010) + decay is referred to as 'slow' due to the small Q-value of the decay. The second selection fully reconstructs the decay D 0 → K 0 S K ± π ∓ , without flavor tagging. In both cases at least one charged particle in the decay chain is required to have a significant impact parameter with respect to any PV.
In the offline selection, trigger signals are associated with reconstructed particles. Selection requirements can therefore be made on the trigger selection itself and on whether the decision was due to the signal candidate, other particles produced in the pp collision, or both. It is required that the hardware hadronic trigger decision is due to the signal candidate, or that the hardware trigger decision is due solely to other particles produced in the pp collision.
Decays K 0 S → π + π − are reconstructed in two different categories: the first involves K 0 S mesons that decay early enough for the pions to be reconstructed in the vertex detector; the second contains K 0 S mesons that decay later such that track segments of the pions cannot be formed in the vertex detector. These categories are referred to as long and downstream, respectively. The long category has better mass, momentum and vertex resolution than the downstream category, and in 2011 was the only category available in the software trigger.
In the simulation, pp collisions are generated using Pythia [21] with a specific LHCb configuration [22]. Decays of hadronic particles are described by EvtGen [23], in which final-state radiation is generated using Photos [24]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [25] as described in Ref. [26].

Signal selection and backgrounds
The offline selection used in this analysis reconstructs the decay chain D * (2010) where the charged pion π + slow from the D * (2010) + decay tags the flavor of the neutral D meson. Candidates are required to pass one of the two software trigger selections described in Sect. 2, as well as several offline requirements. These use information from the RICH detectors to ensure that the charged kaon is well-identified, which reduces the background contribution from the decays D 0 → K 0 S π + π − π 0 and D 0 → K 0 S π − µ + ν µ . In addition the K 0 S decay vertex is required to be well-separated from the D 0 decay vertex in order to suppress the D 0 → K − π + π + π − background, where a π + π − combination is close to the K 0 S mass. D 0 candidates are required to have decay vertices well-separated from any PV, and to be consistent with originating from a PV. This selection suppresses the semileptonic and D 0 → K − π + π + π − backgrounds to negligible levels, while a small contribution from [27] is applied to the reconstructed D * (2010) + decay chain to enhance the resolution in m(K 0 S Kπ), ∆m and the two-body invariant masses m(K 0 S K), m(K 0 S π) and m(Kπ) that are used to probe the resonant structure of these decays. This fit constrains the D * (2010) + decay vertex to coincide with the closest PV with respect to the D * (2010) + candidate, fixes the K 0 S candidate mass to its nominal value, and is required to be of good quality.
Signal yields and estimates of the various background contributions in the signal window are determined using maximum likelihood fits to the m(K 0 S Kπ) and ∆m distributions. The signal window is defined as the region less than 18 MeV/c 2 (0.8 MeV/c 2 ) from the peak value of m(K 0 S Kπ) (∆m), corresponding to approximately three standard deviations of each signal distribution. The three categories of interest are: signal decays, mistagged background where a correctly reconstructed D 0 meson is combined with a charged pion that incorrectly tags the D 0 flavor, and a combinatorial background category, which also includes a small peaking contribution in ∆m from the decay D 0 → K 0 S π + π − π 0 . These fits use candidates in the ranges 139 < ∆m < 153 MeV/c 2 and 1.805 < m(K 0 S Kπ) < 1.925 GeV/c 2 . The sidebands of the m(K 0 S Kπ) distribution are defined as those parts of the fit range where m(K 0 S Kπ) is more than 30 MeV/c 2 from the peak value. The ∆m (m(K 0 S Kπ)) distribution in the signal region of m(K 0 S Kπ) (∆m) is fitted to determine the D * (2010) + (D 0 ) yield in the two-dimensional signal region [28]. The D * (2010) + (D 0 ) signal shape in the ∆m (m(K 0 S Kπ)) distribution is modeled using a Johnson S U [29] (Cruijff [30]) function. In the m(K 0 S Kπ) distribution the combinatorial background is modeled with an exponential function, while in ∆m a power law function is used, f (∆m; m π , p, P, b) = ( ∆m−mπ with the parameters p, P and b determined by a fit in the m(K 0 S Kπ) sidebands. The small D 0 → K 0 S π + π − π 0 contribution in the ∆m distribution is described by a Gaussian function, and the component corresponding to D 0 mesons associated with a random slow pion is the sum of an exponential function and a linear term. These fits are shown in Fig. 2. The results of the fits are used to determine the yields of interest in the two-dimensional signal region. These yields are given in Table 1 for both decay modes, together with the fractions of backgrounds.
A second kinematic fit that also constrains the D 0 mass to its known value is performed and used for all subsequent parts of this analysis. This fit further improves the resolution in the two-body invariant mass coordinates and forces all candidates to lie within the kinematically allowed region of the Dalitz plot. The Dalitz plots [31] for data in the two-dimensional signal region are shown in Fig. 3. Both decays are dominated by a K * (892) ± structure. The K * (892) 0 is also visible as a destructively interfering contribution in the D 0 → K 0 S K − π + mode and the low-m 2 K 0 S π region of the D 0 → K 0 S K + π − mode, while a clear excess is seen in the high-m 2 K 0 S π region. Finally, a veto is applied to candidates close to the kinematic boundaries; this is detailed in Sect. 4.3.

Analysis formalism
The dynamics of a decay D 0 → ABC, where D 0 , A, B and C are all pseudoscalar mesons, can be completely described by two variables, where the conventional choice is to use a pair of squared invariant masses. This paper will use m 2 K 0 S π ≡ m 2 (K 0 S π) and m 2 Kπ ≡ m 2 (Kπ) as this choice highlights the dominant resonant structure of the D 0 → K 0 S K ± π ∓ decay modes.

Isobar models for D
The signal isobar models decompose the decay chain into D 0 → (R → (AB) J )C contributions, where R is a resonance with spin J equal to 0, 1 or 2. Resonances with spin greater than 2 should not contribute significantly to the D 0 → K 0 S K ± π ∓ decays. The corresponding 4-momenta are denoted p D 0 , p A , p B and p C . The reconstructed invariant mass of the resonance is denoted m AB , and the nominal mass m R . The matrix element for the D 0 → K 0 S K ± π ∓ decay is given by Candidates/(100 keV/c 2 ) LHCb Figure 2: Mass (left) and ∆m (right) distributions for the D 0 → K 0 S K − π + (top) and D 0 → K 0 S K + π − (bottom) samples with fit results superimposed. The long-dashed (blue) curve represents the D * (2010) + signal, the dash-dotted (green) curve represents the contribution of real D 0 mesons combined with incorrect π + slow and the dotted (red) curve represents the combined combinatorial and D 0 → K 0 S π + π − π 0 background contribution.  Figure 3: Dalitz plots of the D 0 → K 0 S K − π + (left) and D 0 → K 0 S K + π − (right) candidates in the two-dimensional signal region.
is the momentum of C (A or B) in the R rest frame, and p 0 (q 0 ) is the same quantity calculated using the nominal resonance mass, m R . The meson radius parameters are set to d D 0 = 5.0 (GeV/c) −1 and d R = 1.5 (GeV/c) −1 consistent with the literature [12,33]; the systematic uncertainty due to these choices is discussed in Sect. 5.2. Finally, Ω J (m 2 AB , m 2 AC ) is the spin factor for a resonance with spin J and T R is the dynamical function describing the resonance R. The functional forms for B J (q, q 0 , d) are given in Table 2 and those for Ω J (m 2 AB , m 2 AC ) in Table 3 for J = 0, 1, 2. As the form for Ω 1 is antisymmetric in the indices A and B, it is necessary to define the particle ordering convention used in the analysis; this is done in Table 4. The dynamical function T R chosen depends on the resonance R in question. A relativistic Breit-Wigner form is used unless otherwise noted Several alternative forms are used for specialized cases. The Flatté [34] form is a coupledchannel function used to describe the a 0 (980) ± resonance [12,[35][36][37][38], where the phase space factor is given by and the coupling constants g KK and g ηπ are taken from Ref. [35], fixed in the isobar model fits and tabulated in Appendix A. The Gounaris-Sakurai [39] parameterization is used to describe the ρ(1450) ± and ρ(1700) ± states [37,[40][41][42][43], where and The parameter m K is taken as the mean of m K 0 S and m K ± , and h (m 2 in the limit that m K = m K ± = m K 0 S . Parameters for the ρ resonances ρ(1450) ± and ρ(1700) ± are taken from Ref. [44] and tabulated in Appendix A.
This analysis uses two different parameterizations for the Kπ S-wave contributions, dubbed GLASS and LASS, with different motivations. These forms include both K * 0 (1430) resonance and nonresonant Kπ S-wave contributions. The LASS parameterization takes the form where is an empirical real production form factor, and the phases are defined by The scattering length a, effective range r and K * 0 (1430) mass and width are taken from measurements [45] at the LASS experiment [46] and are tabulated in Appendix A. With the choice f (x) = 1 this form has been used in previous analyses e.g. Refs. [47][48][49], and if δ F is additionally set to zero the relativistic S-wave Breit-Wigner form is recovered. The Watson theorem [50] states that the phase motion, as a function of Kπ invariant mass, is the same in elastic scattering and decay processes, in the absence of final-state interactions (i.e. in the isobar model). Studies of Kπ scattering data indicate that the S-wave remains elastic up to the Kη threshold [45]. The magnitude behavior is not constrained by the Watson theorem, which motivates the inclusion of the form factor f (x), but the LASS parameterization preserves the phase behavior measured in Kπ scattering. The real form factor parameters are allowed to take different values for the neutral and charged K * 0 (1430) resonances, as the production processes are not the same, but the parameters taken from LASS measurements, which specify the phase behavior, are shared between both Kπ channels. A transformed set b = U −1 b of the parameters b = (b 1 , b 2 , b 3 ) are also defined for use in the isobar model fit, which is described in detail in Sect. 5. The constant matrix U is chosen to minimize fit correlations, and the form factor is normalized to unity at the center of the accessible kinematic range, e.g. 1 2 The GLASS (Generalized LASS) parameterization has been used by several recent amplitude analyses, e.g. Refs. [12,37,38], where δ F and δ S are defined as before, and F , φ F and φ S are free parameters in the fit. It should be noted that this functional form can result in phase behavior significantly different to that measured in LASS scattering data when its parameters are allowed to vary freely. This is illustrated in Fig. 10 in Sect. 5.3.

Coherence factor and CP -even fraction
The coherence factor R f and mean strong-phase difference δ f for the multi-body decays D → f and D → f quantify the similarity of the two decay structures [11]. In the limit R f → 1 the matrix elements for the two decays are identical. For D 0 → K 0 S K ± π ∓ the coherence factor and mean strong-phase difference are defined by [10,12] where and the integrals are over the entire available phase space. The restricted phase space coherence factor R K * K e −iδ K * K is defined analogously but with all integrals restricted to an area of phase space close to the K * (892) ± resonance. The restricted area is defined by Ref. [12] as the region where the K 0 S π ± invariant mass is within 100 MeV/c 2 of the K * (892) ± mass. The four observables R K 0 S Kπ , δ K 0 S Kπ , R K * K and δ K * K were measured using quantum-correlated ψ(3770) → D 0 D 0 decays by the CLEO collaboration [12], and the coherence was found to be large for both the full and the restricted regions. This analysis is not sensitive to the overall phase difference between D 0 → K 0 S K + π − and D 0 → K 0 S K + π − . However, since it cancels in δ K 0 S Kπ − δ K * K , this combination, as well as R K 0 S Kπ and R K * K , can be calculated from isobar models and compared to the respective CLEO results.
An associated parameter that it is interesting to consider is the CP -even fraction [51], where the CP eigenstates |D ± are given by 1 √ 2 |D 0 ± |D 0 and B K 0 S Kπ is the ratio of branching fractions of the two D 0 → K 0 S K ± π ∓ modes. As stated above, the relative strong phase δ K 0 S Kπ is not predicted by the amplitude models and requires external input.

Efficiency modeling
The trigger strategy described in Sect. 2, and to a lesser extent the offline selection, includes requirements on variables such as the impact parameter and p T of the various charged particles correlated with the 2-body invariant masses m 2 K 0 S π and m 2 Kπ . There is, Kπ . This efficiency variation is modeled using simulated events generated with a uniform distribution in these variables and propagated through the full LHCb detector simulation, trigger emulation and offline selection. Weights are applied to the simulated events to ensure that various subsamples are present in the correct proportions. These weights correct for known discrepancies between the simulation and real data in the relative reconstruction efficiency for long and downstream tracks, and take into account the ratios of √ s = 7 TeV to √ s = 8 TeV and D 0 → K 0 S K − π + to D 0 → K 0 S K + π − simulated events to improve the description of the data. The efficiencies of offline selection requirements based on information from the RICH detectors are calculated using a data-driven method based on calibration samples [52] of D * (2010) + → D 0 π + slow decays, where D 0 → K − π + . These efficiencies are included as additional weights. A non-parametric kernel estimator [53] is used to produce a smooth function ε(m 2 K 0 S π , m 2 Kπ ) describing the efficiency variation in the isobar model fits. The average model corresponding to the full dataset recorded in 2011 and 2012, which is used unless otherwise noted, is shown in Fig. 4. Candidates very near to the boundary of the allowed kinematic region of the Dalitz plot are excluded, as the kinematics in this region lead to variations in efficiency that are difficult to model. It is required that max(| cos(θ K 0 S π )|, | cos(θ πK )|, | cos(θ KK 0 S )|) < 0.98, where θ AB is the angle between the A and B momenta in the AC rest frame. This criterion removes 5% of the candidates. The simulated events are also used to verify that the resolution in m 2 K 0 S π and m 2 Kπ is around 0.004 GeV 2 /c 4 , corresponding to O(2 MeV/c 2 ) resolution in m(K 0 S π ± ). Although this is not explicitly accounted for in the isobar model fits, it has a small effect which is measurable only on the parameters of the K * (892) ± resonance and is accounted for in the systematic uncertainties.

Fit components
There are three event categories described in Sect. 3 that must be treated separately in the isobar model fits. The signal and mistagged components are described by terms proportional to ε(m 2 , obtained by applying to data in the m(K 0 S Kπ) sidebands the same non-parametric kernel estimator used to model the efficiency variation. The same combinatorial background model is used for both D 0 flavors, and the same efficiency function is used for both modes and D 0 flavors. The overall function used in the fit to D 0 → K 0 S K ± π ∓ decays is therefore Table 1.
All parameters except the complex amplitudes a R e iφ R are shared between the PDFs for both modes and both D 0 flavors. For the other parameters, Gaussian constraints are included unless stated otherwise. The nominal values used in the constraints are tabulated in Appendix A. No constraints are applied for the Kπ S-wave parameters b 1..3 , F , φ S and φ F , as these have no suitable nominal values. The Kπ S-wave parameters a and r are treated differently in the GLASS and LASS models. In the LASS case these parameters are shared between the neutral and charged Kπ channels and a Gaussian constraint to the LASS measurements [45] is included. In the GLASS case these are allowed to vary freely and take different values for the two channels.

Isobar model fits
This section summarizes the procedure by which the amplitude models are constructed, describes the various systematic uncertainties considered for the models and finally discusses the models and the coherence information that can be calculated from them.
Amplitude models are fitted using the isobar formalism and an unbinned maximumlikelihood method, using the GooFit [54] package to exploit massively-parallel Graphics Processing Unit (GPU) architectures. Where χ 2 /bin values are quoted these are simply to indicate the fit quality. Statistical uncertainties on derived quantities, such as the resonance fit fractions, are calculated using a pseudoexperiment method based on the fit covariance matrix.

Model composition
Initially, 15 resonances are considered for inclusion in each of the isobar models: K * (892, 1410, 1680) 0,± , K * 0,2 (1430) 0,± , a 0 (980, 1450) ± , a 2 (1320) ± and ρ(1450, 1700) ± . Pre-liminary studies showed that models containing the K * (1680) resonances tend to include large interference terms, which are cancelled by other large components. Such fine-tuned interference effects are in general unphysical, and are therefore disfavored in the model building [36,55]. The K * (1680) resonances are not considered further, and additionally the absolute value of the sum of interference fractions [56] is required to be less than 30% in all models. In the absence of the K * (1680) resonances, large interference terms are typically generated by the Kπ S-wave contributions. The requirement on the sum of interference fractions, while arbitrary, allows an iterative procedure to be used to search for the best amplitude models. This procedure explores a large number of possible starting configurations and sets of resonances; it begins with the most general models containing all 13 resonances and considers progressively simpler configurations, trying a large number of initial fit configurations for each set of resonances, until no further improvement in fit quality is found among models simple enough to satisfy the interference fraction limit. Higher values of this limit lead to a large number of candidate models with similar fit quality.
A second procedure iteratively removes resonances from the models if they do not significantly improve the fit quality. In this step a resonance must improve the value of −2 log L, where L is the likelihood of the full dataset, by at least 16 units in order to be retained. Up to this point, the K * (892) mass and width parameters and Kπ S-wave parameters have been allowed to vary in the fit, but mass and width parameters for other resonances have been fixed. To improve the quality of fit further, in a third step, S and P -wave resonance parameters are allowed to vary. The tensor resonance parameters are known precisely [3], so remain fixed. At this stage, resonances that no longer significantly improve the fit quality are removed, with the threshold tightened so that each resonance must increase −2 log L by 25 units in order to be retained.
Finally, parameters that are consistent with their nominal values to within 1σ are fixed to the nominal value. The nominal values used are tabulated in Appendix A. The entire procedure is performed in parallel using the GLASS and LASS parameterizations of the Kπ S-wave. The data are found to prefer a solution where the GLASS parameterization of the charged Kπ S-wave has a poorly constrained degree of freedom. The final change to the GLASS models is, therefore, to fix the charged Kπ S-wave F parameter in order to stabilize the uncertainty calculation for the two corresponding a R parameters by reducing the correlations among the free parameters.

Systematic uncertainties
Several sources of systematic uncertainty are considered. Those due to experimental issues are described first, followed by uncertainties related to the amplitude model formalism. Unless otherwise stated, the uncertainty assigned to each parameter using an alternative fit is the absolute difference in its value between the nominal and alternative fit.
As mentioned in Sect. 4.3, candidates extremely close to the edges of the allowed kinematic region of the Dalitz plot are excluded. The requirement made is that the largest of the three | cos(θ AB )| values is less than 0.98. A systematic uncertainty due to this process is estimated by changing the threshold to 0.96, as this excludes a similar additional area of the Dalitz plot as the original requirement.
The systematic uncertainty related to the efficiency model ε(m 2 K 0 S π , m 2 Kπ ) is evaluated in four ways. The first probes the process by which a smooth curve is produced from simulated events; this uncertainty is evaluated using an alternative fit that substitutes the non-parametric estimator with a polynomial parameterization. The second uncertainty is due to the limited sample size of simulated events. This is evaluated by generating several alternative polynomial efficiency models according to the covariance matrix of the polynomial model parameters; the spread in parameter values from this ensemble is assigned as the uncertainty due to the limited sample size. The third contribution is due to possible imperfections in the description of the data by the simulation. This uncertainty is assigned using an alternative simultaneous fit that separates the sample into three categories according to the year in which the data were collected and the type of K 0 S candidate used. As noted in Sect. 2, the sample recorded during 2011 does not include downstream K 0 S candidates. These sub-samples have different kinematic distributions and ε(m 2 K 0 S π , m 2 Kπ ) behavior, so this procedure tests the ability of the simulation process to reproduce the variation seen in the data. The final contribution is due to the reweighting procedures used to include the effect of offline selection requirements based on information from the RICH detectors, and to correct for discrepancies between data and simulation in the reconstruction efficiencies of long and downstream K 0 S candidates. This is evaluated using alternative efficiency models where the relative proportion of the track types is altered, and the weights describing the efficiency of selection requirements using information from the RICH detectors are modified to account for the limited calibration sample size. Additional robustness checks have been performed to probe the description of the efficiency function by the simulated events. In these checks the data are divided into two equally populated bins of the D 0 meson p, p T or η and the amplitude models are re-fitted using each bin separately. The fit results in each pair of bins are found to be compatible within the assigned uncertainties, indicating that the simulated D 0 kinematics adequately match the data.
An uncertainty is assigned due to the description of the hardware trigger efficiency in simulated events. Because the hardware trigger is not only required to fire on the signal decay, it is important that the underlying pp interaction is well described, and a systematic uncertainty is assigned due to possible imperfections. This uncertainty is obtained using an alternative efficiency model generated from simulated events that have been weighted to adjust the fraction where the hardware trigger was fired by the signal candidate.
The uncertainty due to the description of the combinatorial background is evaluated by recomputing the Kπ ) function using m D 0 sideband events to which an alternative kinematic fit has been applied, without a constraint on the D 0 mass. The alternative model is expected to describe the edges of the phase space less accurately, while providing an improved description of peaking features.
An alternative set of models is produced using a threshold of 9 units in the value of −2 log L instead of the thresholds of 16 and 25 used for the model building procedure.
These models contain more resonances, as fewer are removed during the model building process. A systematic uncertainty is assigned using these alternative models for those parameters which are common between the two sets of models.
Two parameters of the Flatté dynamical function, which is used to describe the a 0 (980) ± resonance, are fixed to nominal values in the isobar model fits. Alternative fits are performed, where these parameters are fixed to different values according to their quoted uncertainties, and the largest changes to the fit parameters are assigned as systematic uncertainties.
The effect of resolution in the m 2 K 0 S π and m 2 Kπ coordinates is neglected in the isobar model fits, and this is expected to have an effect on the measured K * (892) ± decay width. An uncertainty is calculated using a pseudoexperiment method, and is found to be small.
The uncertainty due to the yield determination process described in Sect. 3 is measured by changing the fractions f m and f c in the isobar model fit according to their statistical uncertainties, and taking the largest changes with respect to the nominal result as the systematic uncertainty.
There are two sources of systematic uncertainty due to the amplitude model formalism considered. The first is that due to varying the meson radius parameters d D 0 and d R , defined in Sect. 4.1. These are changed from d D 0 = 5.0 (GeV/c) −1 and d R = 1.5 (GeV/c) −1 to 2.5 (GeV/c) −1 and 1.0 (GeV/c) −1 , respectively. The second is due to the dynamical function T R used to describe the ρ(1450, 1700) ± resonances. These resonances are described by the Gounaris-Sakurai functional form in the nominal models, which is replaced with a relativistic P -wave Breit-Wigner function to calculate a systematic uncertainty due to this choice.
The uncertainties described above are added in quadrature to produce the total systematic uncertainty quoted for the various results. For most quantities the dominant systematic uncertainty is due to the meson radius parameters d D 0 and d R . The largest sources of experimental uncertainty relate to the description of the efficiency variation across the Dalitz plot. The fit procedure and statistical uncertainty calculation have been validated using pseudoexperiments and no bias was found.
Tables summarizing the various sources of systematic uncertainty and their relative contributions are included in Appendix C.

Isobar model results
The fit results for the best isobar models using the GLASS and LASS parameterizations of the Kπ S-wave are given in Tables 5 and 6. Distributions of m 2 Kπ , m 2 K 0 S π and m 2 K 0 S K are shown alongside the best model of the D 0 → K 0 S K − π + mode using the GLASS parameterization in Fig. 5. In Fig. 5 and elsewhere the nomenclature R 1 × R 2 denotes interference terms. The corresponding distributions showing the best model using the LASS parameterization are shown in Fig. 6. Distributions for the D 0 → K 0 S K + π − mode are shown in Figs. 7 and 8. Figure 9 shows the GLASS isobar models in two dimensions, and demonstrates that the GLASS and LASS choices of Kπ S-wave parameterization both lead to similar descriptions of the overall phase variation. Figures 5-8 show distributions distorted by efficiency effects, The first uncertainties are statistical and the second systematic.
The first uncertainties are statistical and the second systematic.

LHCb
Total while Fig. 9 shows the decay rate without distortion. Lookup tables for the complex amplitude variation across the Dalitz plot in all four isobar models are available in the supplemental material.
The data are found to favor solutions that have a significant neutral Kπ S-wave contribution, even though the exchange (Fig. 1b) and penguin annihilation (Fig. 1d) processes that contribute to the neutral channel are expected to be suppressed. The expected suppression is observed for the P -wave K * (892) resonances, with the neutral mode

LHCb
Total fit fractions substantially lower. The models using the LASS parameterization additionally show this pattern for the K * (1410) states. The sums of the fit fractions [56], excluding interference terms, in the D 0 → K 0 S K − π + and D 0 → K 0 S K + π − models are, respectively, 103% (109%) and 81% (99%) using the GLASS (LASS) Kπ S-wave parameterization.
Using measurements of the mean strong-phase difference between the D 0 → K 0 S K ± π ∓ modes available from ψ(3770) decays [12], the relative complex amplitudes between each resonance in one D 0 decay mode and its conjugate contribution to the other D 0 decay
Additional information about the models is listed in Appendices B and C, including the interference fractions and decomposition of the systematic uncertainties. The best models also contain contributions from the ρ(1450) ± and ρ(1700) ± resonances in the K 0 S K ± channels, supporting evidence in Ref. [44] of the KK decay modes for these states. Alternative models are fitted where one ρ ± contribution is removed from the best models; in these the value of −2 log L is found to degrade by at least 162 units. Detailed results are tabulated in Appendix B.
The Kπ S-wave systems are poorly understood [6], and there is no clear theoretical guidance regarding the correct description of these systems in an isobar model. As introduced in Sect. 4.1, the LASS parameterization is motivated by the Watson theorem, but this assumes that three-body interactions are negligible and is not, therefore, expected to be precisely obeyed in nature. The isobar models using the GLASS parameterization favor solutions with qualitatively similar phase behavior to those using the LASS parameterization.
S π , m 2 Kπ )| 2 in the best GLASS isobar models, the center row shows the phase behavior of the same models and the bottom row shows the same function subtracted from the phase behavior in the best LASS isobar models. The left column shows the D 0 → K 0 S K − π + mode with D 0 → K 0 S K + π − on the right. The small inhomogeneities that are visible in the bottom row relate to the GLASS and LASS models preferring slightly different values of the K * (892) ± mass and width. Table 7: Modulus and phase of the relative amplitudes between resonances that appear in both the D 0 → K 0 S K − π + and D 0 → K 0 S K + π − modes. Relative phases are calculated using the value of δ K 0 S Kπ measured in ψ(3770) decays [12], and the uncertainty on this value is included in the statistical uncertainty. The first uncertainties are statistical and the second systematic.

Relative amplitude
GLASS LASS 110 ± 20 ± 50 - 70 ± 20 ± 70 -This is illustrated in Fig. 10, which also shows the GLASS forms obtained in fits to D 0 → K 0 S π + π − decays by the BaBar collaboration [37,38] and previously used in fits to the D 0 → K 0 S K ± π ∓ decay modes [12]. This figure shows that the GLASS functional form has substantial freedom to produce different phase behavior to the LASS form, but that this is not strongly favored in the D 0 → K 0 S K ± π ∓ decays. The good quality of fit obtained using the LASS parameterization indicates that large differences in phase behavior with respect The solid (red) curve shows the LASS parameterization, while the dashed (blue) and dash-dotted (green) curves show, respectively, the GLASS functional form fitted to the charged and neutral S-wave channels. The final two curves show the GLASS forms fitted to the charged Kπ S-wave in D 0 → K 0 S π + π − decays in Ref. [37] (triangular markers, purple) and Ref. [38] (dotted curve, black). The latter of these was used in the analysis of D 0 → K 0 S K ± π ∓ decays by the CLEO collaboration [12].
to Kπ scattering data [45,46] are not required in order to describe the D 0 → K 0 S K ± π ∓ decays. A similar conclusion was drawn in Ref.
[57] for the decay D + → K − π + π + , while Ref. [58] found behavior inconsistent with scattering data using the same D + decay mode but a slightly different technique. Ref. [59] studied the Kπ S-wave in τ − → K 0 S π − ν τ decays and found that a parameterization based on the LASS Kπ scattering data, but without a real production form factor, gave a poor description of the τ − decay data.
The quality of fit for each model is quantified by calculating χ 2 using a dynamic binning scheme. The values are summarized in Table 8, while the binning scheme and two-dimensional quality of fit are shown in Appendix B. This binning scheme is generated by iteratively sub-dividing the Dalitz plot to produce new bins of approximately equal population until further sub-division would result in a bin population of fewer than 15 candidates, or a bin dimension smaller than 0.02 GeV 2 /c 4 in m 2 Kπ . This minimum size corresponds to five times the average resolution in these variables.
The overall fit quality is slightly better in the isobar models using the GLASS Kπ S-wave parameterization, but this is not a significant effect and it should be noted that these models contain more degrees of freedom, with 57 parameters fitted in the final GLASS model compared to 50 when using the LASS parameterization. Table 8: Values of χ 2 /bin indicating the fit quality obtained using both Kπ S-wave parameterizations in the two decay modes. The binning scheme for the D 0 → K 0 S K − π + (D 0 → K 0 S K + π − ) mode contains 2191 (2573) bins.

Additional measurements
In this section, several additional results, including those derived from the amplitude models, are presented.

Ratio of branching fractions measurement
The ratio of branching fractions and the restricted region ratio B K * K , defined in the same region near the K * (892) ± resonance as the coherence factor R K * K (Sect. 4.2), are also measured. The efficiency correction due to the reconstruction efficiency ε(m 2 K 0 S π , m 2 Kπ ) is evaluated using the best isobar models, and the difference between the results obtained with the two Kπ S-wave parameterizations is taken as a systematic uncertainty in addition to those effects described in Sect. 5.2. This efficiency correction modifies the ratio of yields quoted in Table 1 by approximately 3%.
The two ratios are measured to be These are the most precise measurements to date.

Coherence factor and CP -even fraction results
The amplitude models are used to calculate the coherence factors R K 0 S Kπ and R K * K , and the strong-phase difference δ K 0 S Kπ − δ K * K , as described in Sect. 4.2. The results are summarized in Table 9, alongside the corresponding values measured in ψ(3770) decays by the CLEO collaboration. Lower, but compatible, coherence is calculated using the isobar models than was measured at CLEO, with the discrepancy larger for the coherence factor calculated over Table 9: Coherence factor observables to which the isobar models are sensitive. The third column summarizes the CLEO results measured in quantum-correlated decays [12], where the uncertainty on δ K 0 S Kπ − δ K * K is calculated assuming maximal correlation between δ K 0 S Kπ and δ K * K .

GLASS LASS CLEO
the full phase space. The results from the GLASS and LASS isobar models are very similar, showing that the coherence variables are not sensitive to the Kπ S-wave parameterization. The coherence factor R K 0 S Kπ , and the ratio of branching fractions B K 0 S Kπ , are combined with the mean phase difference between the two final states measured in ψ(3770) decays [12] to calculate the CP -even fraction F + , defined in Eq. 16, which is determined to be F + = 0.777 ± 0.003 (stat) ± 0.009 (syst), using the GLASS amplitude models. A consistent result is obtained using the alternative (LASS) amplitude models. This model-dependent value is compatible with the direct measurement using only ψ(3770) decay data [12,51].

SU(3) flavor symmetry tests
SU(3) flavor symmetry can be used to relate decay amplitudes in several D meson decays, such that a global fit to many such amplitudes can provide predictions for the neutral and charged K * (892) complex amplitudes in D 0 → K 0 S K ± π ∓ decays [4,5]. Predictions are available for the K * + K − , K * − K + , K * 0 K 0 and K * 0 K 0 complex amplitudes, where K * refers to the K * (892) resonances. There are therefore three relative amplitudes and two relative phases that can be determined from the isobar models, with an additional relative phase accessible if the isobar results are combined with the CLEO measurement of the mean strong phase difference [12]. The results are summarized in Table 10.
The isobar model results are found to follow broadly the patterns predicted by SU(3) flavor symmetry. The amplitude ratio between the K * (892) + and K * (892) 0 resonances, which is derived from the D 0 → K 0 S K − π + isobar model alone, shows good agreement. The two other amplitude ratios additionally depend on the ratio B K 0 S Kπ , and these are more discrepant with the SU(3) predictions. The relative phase between the charged and neutral K * (892) resonances shows better agreement with the flavor symmetry prediction in the D 0 → K 0 S K + π − mode, where both resonances have clear peaks in the data. The GLASS and LASS isobar models are found to agree well, suggesting the problems are not related to the Kπ S-wave.  [5] and results. The uncertainties on phase difference predictions are calculated from the quoted magnitude and phase uncertainties. Note that some theoretical predictions depend on the η-η mixing angle θ η−η and are quoted for two different values. The bottom entry in the table relies on the CLEO measurement [12] of the coherence factor phase δ K 0 S Kπ , and the uncertainty on this phase is included in the statistical uncertainty, while the other entries are calculated directly from the isobar models and relative branching ratio. Where two uncertainties are quoted the first is statistical and the second systematic.

CP violation tests
Searches for time-integrated CP -violating effects in the resonant structure of these decays are performed using the best isobar models. The resonance amplitude and phase parameters a R and φ R are substituted with a R (1 ± ∆a R ) and φ R ± ∆φ R , respectively, where the signs are set by the flavor tag. The convention adopted is that a positive sign produces the D 0 complex amplitude. The full fit results are tabulated in Appendix D.
A subset of the ∆ parameters is used to perform a χ 2 test against the no-CP violation hypothesis: only those parameters corresponding to resonances that are present in the best isobar models using both the GLASS and LASS Kπ S-wave parameterizations are included. The absolute difference |∆ GLASS − ∆ LASS | is assigned as the systematic uncertainty due to dependence on the choice of isobar model. This subset of parameters is shown in Table 11, where the change in fit fraction between the D 0 and D 0 solutions is included for illustrative purposes. In the χ 2 test the statistical and systematic uncertainties are added in quadrature.

Conclusions
The decay modes D 0 → K 0 S K ± π ∓ have been studied using unbinned, time-integrated, fits to a high purity sample of 189 670 candidates, and two amplitude models have been constructed for each decay mode. These models are compared to data in a large number of bins in the relevant Dalitz plots and a χ 2 test indicates a good description of the data.
Models are presented using two different parameterizations of the Kπ S-wave systems, which have been found to be an important component of these decays. These systems are poorly understood, and comparisons have been made to previous results and alternative parameterizations, but the treatment of the Kπ S-wave is found to have little impact on the other results presented in this paper. The large fractions attributed to the neutral Kπ S-wave channels could indicate larger than expected contributions from the penguin annihilation diagrams shown in Fig. 1d.
The models are seen to favor small, but significant, contributions from the ρ(1450, 1700) ± → K 0 S K ± resonances, modes which were seen by the OBELIX experiment [44] but are not well established. All models contain clear contributions from both the K * (892) ± and K * (892) 0 resonances, with the K * (892) 0 contribution found to be suppressed as expected from the diagrams shown in Fig. 1. This allows the full set of amplitudes in these decays that are predicted by SU(3) flavor symmetry to be tested, in contrast to the previous analysis by the CLEO collaboration [12]. Partial agreement is found with these predictions.
The ratio of branching fractions between the two D 0 → K 0 S K ± π ∓ modes is also measured, both across the full Dalitz plot area and in a restricted region near the K * (892) ± resonance, with much improved precision compared to previous results.
Values for the D 0 → K 0 S Kπ coherence factor are computed using the amplitude models, again both for the whole Dalitz plot area and in the restricted region, and are found to be in reasonable agreement with direct measurements by CLEO [12] using quantum-correlated ψ(3770) → D 0 D 0 decays. The CP -even fraction of the D 0 → K 0 S Kπ decays is also computed, using input from the quantum-correlated decays, and is found to be in agreement with the direct measurement [12,51]. A search for time-integrated CP violation is carried out using the amplitude models, but no evidence is found with either choice of parameterization for the Kπ S-wave.
The models presented here will be useful for future D 0 -D 0 mixing, indirect CP violation and CKM angle γ studies, where knowledge of the strong-phase variation across the Dalitz plot can improve the attainable precision. These improvements will be particularly valuable for studies of the large dataset that is expected to be accumulated in Run 2 of the LHC.

Appendices A Additional isobar formalism information
This appendix contains Table 12, which summarizes the nominal values used for the various resonance and form factor parameters.

B Additional isobar model information
This appendix contains additional information about the various isobar model parameters that are used and allowed to vary freely in the model fits, e.g. resonance mass and width parameter values, and parameters of the GLASS and LASS Kπ S-wave functional forms.
Tables 13 and 14 summarize the most significant interference terms in the D 0 → K 0 S K − π + and D 0 → K 0 S K + π − models, respectively. Table 15 defines the matrices U used to define the LASS Kπ S-wave form factor. Tables 16 (GLASS) and 17 (LASS) summarize the various resonance and form factor parameters. The nominal values that are used in Gaussian constraint terms are given in Appendix A. Figure 11 shows the smooth functions that describe the combinatorial background in the isobar model fits. Figure 12 illustrates the two-dimensional quality of fit achieved in the four isobar models and shows the binning scheme used to derive χ 2 /bin values.
The changes in −2 log L obtained in alternative models where one ρ contribution is removed are given in Table 18.
S π ), used to describe the combinatorial background component in the D 0 → K 0 S K − π + (left) and D 0 → K 0 S K + π − (right) amplitude model fits.

Parameter
Value

Parameter
Value Parameter Value MeV/c 2 ρ(1700) ± m R 1552 ± 13 ± 26 MeV/c 2 Table 18: Change in −2 log L value when removing a ρ resonance from one of the models.

C Systematic uncertainty tables
This appendix includes tables summarizing the various contributions to the systematic uncertainties assigned to the various results. The table headings correspond to the uncertainties discussed in Sect. 5.2 with some abbreviations to allow the tables to be typeset compactly. Definitions of the various abbreviations are given in Table 19. The quantity 'DFF' listed in the tables is the sum of fit fractions from the various resonances, excluding interference terms. Tables 20 (GLASS) and 21 (LASS) show the results for the complex amplitudes and fit fractions in the D 0 → K 0 S K − π + models, Tables 22 (GLASS) and 23 (LASS) show the corresponding values for the D 0 → K 0 S K + π − models and Tables 24 and 25 summarize the uncertainties for the parameters that are not specific to a decay mode.
In each of these tables the parameter in question is listed on the left, followed by the central value and the corresponding statistical (first) and systematic (second) uncertainty. The subsequent columns list the contributions to this systematic uncertainty, and are approximately ordered in decreasing order of significance from left to right.

D CP violation fit results
This appendix contains Table 26, which summarizes the full fit results of the CP violation searches described in Sect. 6.4.

Abbreviation Description max(| cos |)
Variation of the cut that excludes the boundary regions of the Dalitz plot. Efficiency Two efficiency modelling uncertainties added in quadrature: using an alternative parameterization, and accounting for the limited size of the simulated event sample. Joint Uncertainty obtained by simultaneously fitting disjoint sub-sets of the dataset, separated by the year of data-taking and type of K 0 S daughter track, with distinct efficiency models. Weights Three uncertainties related to the re-weighting of simulated events used to generate the efficiency model ε(m 2 K 0 S π , m 2 Kπ ), added in quadrature. These account for: incorrect simulation of the underlying pp interaction, uncertainty in the relative yield of long and downstream K 0 S candidates, and uncertainty in the efficiency of selection requirements using information from the RICH detectors. Comb.
Using an alternative combinatorial background model. −2 log L Using a more complex alternative model where the threshold in ∆(−2 log L) for a resonance to be retained is reduced to 9 units. Flatté Variation of the Flatté lineshape parameters for the a 0 (980) ± resonance according to their nominal uncertainties. f m , f c Variation of the mistag and combinatorial background rates according to their uncertainties in the mass fit.
Variation of the meson radius parameters. T ρ ± Switching to a Breit-Wigner dynamical function to describe the ρ(1450, 1700) ± resonances.
zation. The headings are defined in Table 19.
zation. The headings are defined in Table 19.
zation. The headings are defined in Table 19.

Var
Baseline  Table 24: Systematic uncertainties for shared parameters, coherence and relative branching ratio observables in the GLASS models.
The systematic uncertainty on the K * (892) ± width due to neglecting resolution effects in the nominal models is 0.5 MeV/c 2 MeV/c 2 MeV/c 2 Table 25: Systematic uncertainties for shared parameters, coherence and relative branching ratio observables in the LASS models.
The systematic uncertainty on the K *

Supplemental material
This is divided into two parts: lookup tables for the complex amplitude and covariance information, each for the four quoted amplitude models. These are available at Ref. [60]. No correlation information is included for systematic uncertainties.

Isobar model lookup tables
The lookup table filenames are listed in Table 27. As an example, the first five lines of the file glass fav lookup.txt are: The first three lines are comments, describing which D 0 decay mode and isobar model this file corresponds to, giving the precise nominal masses used in the fit and, finally, defining the data fields in the remainder of the file. As this shows, the models are evaluated on a grid with a spacing of 0.00125 GeV 2 /c 4 .

Covariance information
A reduced covariance matrix is presented for each isobar model, tabulating the correlations between the complex amplitudes a R e iφ R . These are listed in files named analogously to those in Table 27, e.g. glass fav covariance.txt. An example first four lines: