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

Amplitude models are constructed to describe the resonance structure of ${D^{0}\to K^{-}\pi^{+}\pi^{+}\pi^{-}}$ and ${D^{0} \to K^{+}\pi^{-}\pi^{-}\pi^{+}}$ decays using $pp$ collision data collected at centre-of-mass energies of 7 and 8 TeV with the LHCb experiment, corresponding to an integrated luminosity of $3.0\mathrm{fb}^{-1}$. The largest contributions to both decay amplitudes are found to come from axial resonances, with decay modes $D^{0} \to a_1(1260)^{+} K^{-}$ and $D^{0} \to K_1(1270/1400)^{+} \pi^{-}$ being prominent in ${D^{0}\to K^{-}\pi^{+}\pi^{+}\pi^{-}}$ and $D^{0}\to K^{+}\pi^{-}\pi^{-}\pi^{+}$, respectively. Precise measurements of the lineshape parameters and couplings of the $a_1(1260)^{+}$, $K_1(1270)^{-}$ and $K(1460)^{-}$ resonances are made, and a quasi model-independent study of the $K(1460)^{-}$ resonance is performed. The coherence factor of the decays is calculated from the amplitude models to be $R_{K3\pi} = 0.459\pm 0.010\,(\mathrm{stat}) \pm 0.012\,(\mathrm{syst}) \pm 0.020\,(\mathrm{model})$, which is consistent with direct measurements. These models will be useful in future measurements of the unitary-triangle angle $\gamma$ and studies of charm mixing and $C\!P$ violation.


Introduction
The decays 1 D 0 → K − π + π + π − and D 0 → K + π − π + π − have an important role to play in improving knowledge of the Cabibbo-Kobayashi-Maskawa (CKM) unitarity-triangle angle γ ≡ arg(−V ud V * ub /V cd V * cb ). Sensitivity to this parameter can be obtained by measuring CP -violating and associated observables in the decay B − → DK − , where D indicates a neutral charm meson reconstructed in final states common to both D 0 and D 0 , of which K ∓ π ± π ± π ∓ are significant examples [1,2]. A straightforward approach to such an analysis is to reconstruct the four-body D-meson decays inclusively, which was performed by the LHCb collaboration in a recent measurement [3]. Alternatively, additional sensitivity can be sought by studying the variation of the observables across the phase space of the D decays, a strategy that requires knowledge of the variation of the decay amplitudes of the charm mesons.
Studies of charm mixing and searches for CP violation in the D 0 -D 0 system, which for these final states have only been performed inclusively [4], will also benefit from an understanding of the variation of the decay amplitudes across their phase space. These decay modes are also a rich laboratory for examining the behaviour of the strong interaction at low energy, through studies of the intermediate resonances that contribute to the final states. All these considerations motivate an amplitude analysis of the two decays.
The decay D 0 → K − π + π + π − has a branching ratio of (8.29 ± 0.20)% [5], which is the highest of all D 0 decay modes involving only charged particles, and is predominantly mediated by Cabibbo-favoured (CF) transitions. The decay D 0 → K + π − π − π + is dominated by doubly Cabibbo-suppressed (DCS) amplitudes, with small contributions from mixingrelated effects, and occurs at a rate that is suppressed by a factor of (3.22 ± 0.05) × 10 −3 [4] compared to that of the favoured mode. The favoured and suppressed modes are here termed the 'right-sign' (RS) and 'wrong-sign' (WS) decay, respectively, on account of the charge correlation between the kaon and the particle used to tag the flavour of the parent meson.
In this paper, time-integrated amplitude models of both decay modes are constructed using pp collision data collected at centre-of-mass energies of 7 and 8 TeV with the LHCb experiment, corresponding to an integrated luminosity of 3.0 fb −1 . The RS sample size is around 700 times larger than the data set used by the Mark III collaboration to develop the first amplitude model of this decay [6]. An amplitude analysis has also been performed on the RS decay by the BES III collaboration [7] with around 1.6% of the sample size used in this analysis. This paper reports the first amplitude analysis of the WS decay.
The paper is organised as follows. In Sect. 2 the detector, data and simulation samples are described, and in Sect. 3 the signal selection is discussed. The amplitude-model formalism is presented in Sect. 4, and the fit method and model-building procedure in Sect. 5. Section 6 contains the fit results and conclusions are drawn in Sect. 7.

Detector and simulation
The LHCb detector [8] 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 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. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
The trigger [9] 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. The software trigger requires a two-, three-or four-track secondary vertex with a significant displacement from the primary pp interaction vertices. At least one charged particle must have a transverse momentum p T > 1.7(1.6) GeV/c and be inconsistent with originating from a PV. A multivariate algorithm [10] is used for the identification of secondary vertices consistent with the decay of a b hadron.
In the simulation, pp collisions are generated using Pythia [11] with a specific LHCb configuration [12]. Particle decays are described by EvtGen [13]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [14] as described in Ref. [15].

Signal selection and backgrounds
The decay chain B → D * (2010) + µ − X with D * (2010) + → D 0 π + slow is reconstructed as a clean source of D 0 mesons for analysis. The D 0 mesons are reconstructed in the K ∓ π ± π ± π ∓ final states. The charged pion, π + slow , originating from the D * (2010) + is referred to as 'slow' due to the small Q-value of the decay. The charge of the muon and slow pion are used to infer the flavour of the neutral D meson. Candidates are only accepted if these charges lead to a consistent hypothesis for the flavour of the neutral D meson. All other aspects of the reconstruction and selection criteria are identical between the RS and WS samples.
The two-dimensional plane m Kπππ vs. ∆m, where m Kπππ is the invariant mass of the D 0 meson candidate, and ∆m = m Kππππ slow −m Kπππ is mass difference between the D * (2010) + and D 0 meson candidates, is used to define signal and sideband regions with which to perform the amplitude analysis and study sources of background contamination. The signal region is defined as ±0.75 MeV/c 2 (±18 MeV/c 2 ) of the signal peak in ∆m(m Kπππ ), which corresponds to about three times the width of the peak.
It is required that the hardware trigger decision is either due to the muon candidate or is independent of the particles constituting the reconstructed decay products of the B candidate. For example, a high-p T particle from the other B meson decay in the event firing the hadron trigger. The software trigger decision is required to either be due to the muon candidate or a two-three-or four-track secondary secondary vertex.
The WS sample is contaminated by a category of RS decays in which the kaon is mis-identified as a pion, and a pion as a kaon. To suppress this background, it is required that the kaon is well identified by the RICH detectors. The residual contamination from this background is removed by recalculating the mass of the D 0 candidate with the mass hypotheses of a kaon and each oppositely charged pion swapped, then vetoing candidates that fall within 12 MeV/c 2 of the nominal mass of the D 0 meson. As the majority of particles from the PV are pions, the particle identification requirements on the kaon also reduces the background from random combinations of particles.
Remaining background from random combinations of particles can be divided into two categories. Candidates where the D 0 is reconstructed from a random combination of tracks are referred to as combinatorial background. Candidates where the D 0 is correctly reconstructed but paired with an unrelated π + slow are referred to as mistag background. This latter source of background is dominated by RS decays. Both of these backgrounds are suppressed using a multivariate classifier based on a Boosted Decision Tree (BDT) [16][17][18] algorithm. The BDT is trained on RS data candidates from the signal region and the sidebands of the WS data, and uses 15 variables related to the quality of the reconstruction of the PV, B and D 0 decay vertices, and the consistency of tracks in the signal candidate incoming from these vertices. Variables pertaining to the D 0 kinematics and its decay products are avoided to minimise any bias of the phase-space acceptance.
The signal and background yields in the signal region for each sample is determined by simultaneously fitting the two-dimensional ∆m vs. m Kπππ distribution for both samples. The D 0 , muon and slow pion candidates are constrained to originate from a common vertex in calculating the D 0 and D * + masses. This requirement improves the resolution of the ∆m distribution by approximately a factor of two. The signal is modelled with a product of two Cruijff [19] functions. The Cruijff shape parameters are shared between both samples. The combinatorial background is modelled by a first-order polynomial in m Kπππ , and by a threshold function in ∆m, where Q = ∆m − m π and the parameters p, a are determined by the fit. The background shape parameters, including those for the polynomial in m Kπππ , are allowed to differ between WS and RS samples. The mistag background component is a product of the signal shape in m Kπππ and the combinatorial background shape in ∆m. The optimal requirement on the output of the BDT classifier is selected by repeating the fit varying this requirement, and maximising the expected significance of the WS signal, which is defined as where N bkg is the background yield in the signal region. The expected number of WS candidates,N sig , is estimated by scaling the number of RS signal candidates in the signal region by the ratio of branching fractions. The yields of the various contributions for both samples are listed in Table 1, and the m Kπππ and ∆m distributions, with the fit projections superimposed, are shown in Fig. 1   arising from mistagged decays. Studies of simulated data indicate that the selected sample has a relatively uniform acceptance across the phase space, with approximately 30% reductions in acceptance near the edges of the kinematically allowed region. The samples also have a relatively uniform selection efficiency in decay time, being constant within ±10% for lifetimes greater than one average lifetime of the D meson.
For the amplitude analysis, a kinematic fit is performed constraining the D 0 mass to its known value [20], which improves the resolution in the D 0 phase space. This also forces all candidates to lie inside the kinematically allowed region. Candidates are only accepted if this kinematic fit converges.

Formalism of amplitude model
The amplitudes contributing to the decays D 0 → K ∓ π ± π ± π ∓ are described in terms of a sequence of two-body states. It is assumed that once these two-body states are produced, rescattering against other particles can be neglected. Two-body processes are  [21,22] is adopted, and will be denoted by [π + π − ] L=0 and [K ∓ π ± ] L=0 throughout for π + π − and K ∓ π ± S-waves, respectively.
The following decay chains are considered: Cascade decays have the topology D 0 → X [Y [P 1 P 2 ] P 3 ] P 4 -the D 0 meson decays into a stable pseudoscalar state P 4 and an unstable state X. The unstable state then decays to three pseudoscalars P 1,2,3 via another intermediate unstable state (Y ). There are three distinct possibilities for cascade decays. The resonance X can either have isospin I = 1/2, and will therefore decay into the K ∓ π ± π ∓ final state, or have isospin I = 1 and therefore will decay into the π + π − π ± final state. In the K ∓ π ± π ∓ case, the next state in the cascade Y can either be in K ∓ π ± or π + π − , referred to as cases (1) and (2), respectively. In the π + π − π ± case, there is only the π + π − state, referred to as case (3).
Two complex parameters can be used to describe cascade decays: the coupling between the D 0 meson and the first isobar, and then the coupling between the first isobar and the second intermediate state. One of the couplings between isobars can be fixed by convention, typically the dominant channel. For example, for the a 1 (1260) + resonance, the couplings for subdominant decay chains such as a 1 (1260) + → [π + π − ] L=0 π + are defined with respect to the dominant a 1 (1260) + → ρ(770) 0 π + decay.
Quasi two-body decays have the topology D 0 → X [P 1 P 2 ] Y [P 3 P 4 ] -the D 0 meson decays into a pair of unstable states, which in turn each decay to a pair of stable pseudoscalar mesons. The only possibility where X, Y form resonances of conventional quark content is X [K − π + ] Y [π + π − ], with an example of a typical process being The parameters to be determined describe the coupling between the D 0 initial state and the quasi two-body state. In the above example, there are three different possible orbital configurations of the vector-vector system, and hence this component has three complex parameters.
Decay chains are described using a product of dynamical functions for each isobar and a spin factor. The amplitude for each decay chain is explicitly made to respect Bose symmetry by summing over both possible permutations of same-sign pions. The total amplitude is then modelled as a coherent sum of these processes. Spin factors are modelled using the Rarita-Schwinger formalism following the prescription in Ref. [23]; the details of this formulation are included in Appendix A.
Resonances are modelled with the relativistic Breit-Wigner function unless otherwise stated, which as a function of the invariant-mass squared, s, takes the form where the mass of the resonance is m 0 and Γ(s) is the energy-dependent width. The form factor for a decay in which the two decay products have relative orbital angular momentum L is given by the normalised Blatt-Weisskopf function [24] B L (q, 0), where q is the three-momentum of either decay product in the rest frame of the resonance, and is normalised to unity at zero momentum transfer. The factor k normalises the lineshape integrated over all values of s if the Blatt-Weisskopf form-factor and energy dependence of the width are neglected, and is included to reduce correlations between the coupling to the channel and the mass and width of the resonance. For a resonance that decays via a single channel to two stable particles, such as ρ(770) 0 → π + π − , the width is given by where Γ 0 is the width at the resonance mass, and q 0 is the linear momentum of either decay product evaluated at the rest mass of the resonance. The energy-dependent width of a resonance that decays to a three-body final-state must account for the dynamics of the intermediate decay process, and follows that developed for the decay τ + → a 1 (1260) + ν τ by the CLEO collaboration in Ref. [25]. The width of a resonance R decaying into three bodies abc can be expressed in terms of the spin-averaged matrix element of the decay M R→abc integrated over the phase space of the three-particle final state, where the matrix element consists of a coherent sum over the intermediate states in the three-body system, described using the isobar model and using the fitted couplings between the resonance and the intermediate isobars.
In the example of the decay of the a 1 (1260) + resonance, these are predominately the couplings to the ρ(770) 0 π + and [π + π − ] L=0 π + intermediate states. The width is normalised such that Γ(m 2 0 ) = Γ 0 . In the three-body case, exponential form-factors are used rather than normalised Blatt-Weisskopf functions, where r characterises the radius of the decaying resonance. The K-matrix formalism [21] provides a convenient description of a two-particle scattering amplitude, which is particularly useful in parameterising S-wave systems. This formulation can then be used in the description of multibody decays on the assumption that rescattering against the other particles in the decay can be neglected. The K-matrix formalism is used in this analysis to describe the π + π − and K ∓ π ± S-waves due to its relative success in parameterising the scalar contributions to three-body decays [26,27] of the D meson.
The π + π − S-wave (isoscalar) amplitude is modelled using the K matrix from Ref. [26,28], which describes the amplitude in the mass range 280 MeV/c 2 < √ s < 1900 MeV/c 2 , considering the effects of five coupled channels, ππ, KK, ππππ, ηη, η η, and five poles with masses which generate the resonances. The K matrix also includes polynomial terms that describe nonresonant scattering between hadrons. The coupling to each of these poles and the direct coupling to each of the five channels depend on the production mode, which is modelled using the production vector or P-vector approach, in which the amplitude is whereρ is the two-body phase-space matrix. The complex-valued vector function,P , has one component for each of the coupled channels, and describes the coupling between the initial state and either one of the poles or a direct coupling to one of these channels. The generic P-vector for the isoscalar K-matrix therefore has 10 complex parameters. An additional complexity in the four-body case is that there are several initial states that couple to the π + π − S-wave, each of which has its own P vector. Several simplifying assumptions are therefore made to the P vector to avoid introducing an unreasonable number of degrees of freedom. The only direct production terms included in the P vector are to the ππ and KK states, as the production of the π + π − final state via a direct coupling to another channel all have similar structure below their corresponding production thresholds. The couplings to poles 3, 4 and 5 (where the numbering of the poles is defined in Ref. [28]) are also fixed to zero, as production of these poles only has a small effect within the phase space. This choice reduces the number of free parameters per S-wave production mechanism to four complex numbers. The couplings to the poles are described by β 0 and β 1 , while the direct couplings to each channel by f ππ and f KK .
The production vectors used here should therefore be considered as a minimal simplified model. For production of π + π − S-wave states via resonances, such as the decay chain a 1 (1260) + → [π + π − ] L=0 π + , improved sensitivity to the structure of the π + π − state can be achieved by studying a decay mode that produces the a 1 (1260) + with a larger phase space. In several cases, one or more of these couplings are found to be negligible for a given production mode, and therefore are fixed to zero.
The K ∓ π ± S-wave is modelled using the K matrices from the analysis of D + → K − π + π + by the FOCUS collaboration [27]. The I = 1/2 K matrix considers two channels, Kπ and Kη , and a single pole which is responsible for generating the K * (1430) 0 resonance. Additionally, the K matrix includes polynomial terms that describe nonresonant scattering between the hadrons. The K ∓ π ± S-wave also contains a I = 3/2 component. No poles or inelasticity are expected with this isospin, and therefore the associated amplitude can be modelled using a K matrix consisting of a single scalar term.
The I = 1/2 amplitudes are constructed in the Q-vector [22] approximation. The P vector has the same pole structure as the K matrix, and therefore the approximation can be made, whereα(s) is a slowly varying complex vector. This is sometimes referred to as the Q-vector [22] approximation, and allows the insertion ofK −1K into Eq. 7, and the rephrasing of the I = 1/2 decay amplitude, A 1/2 , in terms of the T-matrix elements from scattering: which is the transition matrix associated with the I = 1/2 scattering process. Given the relatively small energy range available to the K ∓ π ± system, it is reasonable to approximatê α(s) as a constant. Inclusion of polynomial terms inα(s) is found not to improve the fit quality significantly. The coupling to the Kη channel, α Kη , is defined with respect to the coupling to the Kπ channel, α Kπ in all production modes. If the phase of α Kη is zero, the phase shift of the I = 1/2 component matches that found in scattering experiments, which is the expected result if Watson's theorem [29] holds for these decays. Similar to the π + π − S-wave, the components ofα and the coupling to the I = 3/2 channel are allowed to differ between production modes.

Fit formalism and model construction
Independent fits are performed on the D 0 → K − π + π + π − and D 0 → K + π − π − π + data sets, using an unbinned maximum likelihood procedure to determine the amplitude parameters. The formalism of the fit is described in Sects. 5.1-5.3, and the method for systematically selecting plausible models is discussed in Sect. 5.4.

Likelihood
The probability density functions (PDFs) are functions of position in D 0 decay phase-space, x, and are composed of the signal amplitude model and the two sources of background described in Sect. 3: The signal PDF is described by the function |M(x)| 2 , where M(x) is the total matrix element for the process, weighted by the four-body phase-space density φ(x), and the phase-space acceptance, ε(x). The mistag component involving M(x), is only present in the WS sample, and is modelled using the RS signal PDF. The combinatorial background is modelled by P c (x), and is present in both samples. The normalisation of each component is given by the integral of the PDF over the phase space, N i , where i = (c, s, m), weighted by the fractional yield, Y i , determined in Sect. 3. The PDF that describes the combinatorial background in the WS sample is fixed to the results of a fit to the two sidebands of the m Kπππ distribution, below 1844.5 MeV/c 2 and above 1888.5 MeV/c 2 . The components in this model are selected using the same algorithm to determine the resonant content of the signal modes, which is discussed in Sect. 5.4. In this case, the PDF incoherently sums the different contributions and assumes no angular correlations between tracks. The contamination from combinatorial background in the RS sample is very low, and hence this contribution can safely be assumed to be distributed according to phase space, that is P c (x) = 1.
The function to minimise is As the efficiency variation across the phase space factorises in the PDF, these variations result in a constant shift in the likelihood everywhere except the normalisation integrals, and hence can be neglected in the minimisation procedure. Efficiency variations can then be included in the fit by performing all integrals using simulated events that have been propagated through the full LHCb detector simulation and selection. These events are referred to as the integration sample. The values of the normalisation integrals are independent of the generator distribution of the integration sample, however the uncertainties on the integrals are minimised when integration events approximate the function being integrated, which is known as importance sampling. Therefore, integration samples are generated using preliminary models that do not include efficiency effects.

Goodness of fit
The quality of fits is quantified by computing a χ 2 metric. Candidates are binned using an adaptive binning scheme. Five coordinates are selected, and the phase space is repeatedly divided in these coordinates such that each bin contains the same number of candidates, following the procedure described in Ref. [4]. The division is halted when each bin contains between 10 and 20 entries. This procedure results in 32,768 approximately equally populated bins for the RS sample, and 256 for the WS sample. Five two-and three-body invariant mass-squared combinations are used as coordinates for the binning procedure, s π + π − π + , s K − π + , s K − π − , s π + π − and s K − π + π − . The χ 2 is defined as where N i is the observed number of candidates in bin i and N i , the expected number of entries determined by reweighting the integration sample with the fitted PDF. The statistical uncertainty from the limited size of the integration sample,σ i , is included in the definition of the χ 2 , and is estimated as where ω j is the weight of integration event j. The χ 2 per degree of freedom is used as the metric to optimise the decay chains included in a model, using the model-building procedure described in Sect. 5.4.

Fit fractions
The values of coupling parameters depend strongly on various choices of convention in the formalism. Therefore, it is common to define the fractions in the data sample associated with each component of the amplitudes (fit fractions). In the limit of narrow resonances, the fit fractions are analogous to relative branching fractions. The fit fraction for component i is For cascade processes, the different secondary isobars contribute coherently to the fit fractions. The partial fit fractions for each sub-process are then defined as the fit fraction with only the contributions from the parent isobar included in the denominator.

Model construction
The number of possible models that could be used to fit the amplitudes is extremely large due to the large number of possible decay chains (≈ 100). A full list of the components considered is included in Appendix B. A model of "reasonable" complexity typically contains O(10) different decay chains. Therefore, the number of possible models is extremely large, and only an infinitesimal fraction of these models can be tested. An algorithmic approach to model building is adopted, which begins with an initial model and attempts to iteratively improve the description by adding decay chains. For D 0 → K − π + π + π − the initial model is that constructed by the Mark III collaboration [6], augmented by knowledge from other analyses, such as the additional decay channels of the a 1 (1260) + found in the amplitude analysis of the decay D 0 → π + π − π + π − performed by the FOCUS collaboration [30]. The two-body nonresonant terms in the Mark III model are replaced with the relevant K matrices, and the four-body nonresonant term replaced with a quasi two-body scalar-scalar term [K − π + ] L=0 [π + π − ] L=0 , modelled using a product of K matrix amplitudes.
For D 0 → K + π − π − π + , where no previous study exists, the initial model is obtained by inspection of the invariant-mass distributions. There are clear contributions from the K * (892) 0 and ρ(770) 0 resonances, and therefore combined with the expectation that the vector-vector contributions should be similar between WS and RS, the quasi two-body mode D 0 → K * (892) 0 ρ(770) 0 is included in all three allowed orbital states L = (0, 1, 2). The scalar-scalar contribution should also be comparable between WS and RS decay modes, and hence the quasi two-body term The steps of the model-building procedure are 1. Take a model and a set of possible additional decay chains, initially the complete set discussed in Appendix. B. Perform a fit to the data using this model adding one of these decay chains.
2. If adding this decay chain improves the χ 2 per degree of freedom by at least 0.02, then retain the model for further consideration.
3. On the first iteration, restrict the pool of decay chains that are added to the model to those 40 contributions that give the largest improvements to the fit.
4. Reiterate the model-building procedure, using the 15 models with the best fit quality from step 2 as starting points. Finish the procedure if no model has improved significantly.
The model-building procedure therefore results in an ensemble of parametrisations of comparable fit quality.

Fit results
This section presents fit results and systematic uncertainties, with the latter discussed first in Sect. 6.1. The model-building procedure outlined in Sect. 5.4 results in ensembles of parameterisations of comparable fit quality. The models discussed in this section, which are referred to as the baseline models, and are built to include all decay chains that are common to the majority of models that have a χ 2 per degree of freedom differing from the best-fitting models by less than 0.1. The results for these baseline models are shown and their features discussed in Sect. 6.2 and Sect. 6.3 for the RS decay and the WS decay, respectively. The general features of models in the ensembles are discussed in Sect. 6.4. In Sect. 6.5 the models are used to calculate the coherence factor of the decays, and an assessment is made of the stability of the predicted coherence factors, strong-phase differences and amplitude ratios with respect to the choice of WS model in regions of phase space.

Systematic uncertainties
Several sources of systematic uncertainty are considered. Experimental issues are discussed first, followed by uncertainties related to the model and the formalism. All parameters in the fit have a systematic uncertainty originating from the limited size of the integration sample used in the likelihood minimisation. This effect is reduced by importance sampling. The remaining uncertainty is estimated using a resampling technique. Half of the integration sample is randomly selected, and the fit performed using only this subsample. This procedure is repeated many times, and the systematic uncertainty from the finite integration statistics is taken to be 1/ √ 2 of the spread in fit parameters.
There is an additional systematic uncertainty due to the imperfect simulation, which affects the efficiency corrections. The RS data are divided into bins in the D 0 transverse momentum, in which the efficiency corrections may be expected to vary, and the fit is performed independently in each bin. The results of these fits are combined in an uncertainty-weighted average, including the correlations between the different parameters, and the absolute difference between the parameters measured by this procedure and the usual fitting procedure is assigned as the systematic uncertainty. Additionally, the data is divided by data-taking year and software trigger category and independent fits performed using these subsamples. The fit results are found to be compatible within the assigned uncertainties between these samples, hence no additional systematic uncertainty is assigned.
The uncertainty associated with the determination of the signal fraction and mistag fraction in each sample is measured by varying these fractions within the uncertainties found in the fit to the m Kπππ vs. ∆m plane.
Parameters that are fixed in the fit, such as the ρ(770) 0 mass and width, are randomly varied according to the uncertainties given in Ref. [20], and the corresponding spreads in fit results are assigned as the uncertainties. It is assumed that input correlations between these parameters are negligible. When performing fits to the WS sample, several parameters, such as the mass, width and couplings of the K 1 (1270) ± resonance, are fixed to the values found in the RS fit. The uncertainty on these parameters is propagated to the WS fit by randomly varying these parameters by their uncertainties. The radii of several particles used in the Blatt-Weisskopf form factor are varied using the same procedure. The D 0 radial parameter is varied by ±0.5 GeV −1 c.
The uncertainty due to the background model in the WS fit is estimated using pseudoexperiments. A combination of simulated signal events generated with the final model and candidates from outside of the D 0 signal region is used to approximate the real data. The composite dataset is then fitted using the signal model, and differences between the true and fitted values are taken as the systematic uncertainties on the background parametrisation.
The choice of model is an additional source of systematic uncertainty. It is not meaningful to compare the coupling parameters between different parametrisations, as these are by definition the parameters of a given model. It is however useful to consider the impact the choice of parametrisation has on fit fractions and the fitted masses and widths. Therefore, the model choice is not included in the total systematic uncertainty, but considered separately in Sect. 6.4-6.5.
The total systematic uncertainty is obtained by summing the components in quadrature. The total systematic uncertainty is significantly larger than the statistical uncertainty on the RS fit, with the largest contributions coming from the form factors that account for the finite size of the decaying mesons. For the WS fit, the total systematic uncertainty is comparable to the statistical uncertainty, with the largest contribution coming from the parametrisation of the combinatorial background. A full breakdown of the different sources of systematic uncertainty for all parameters is given in Appendix C.

Results for the RS decay
Invariant mass-squared projections for D 0 → K − π + π + π − are shown in Fig. 2 together with the expected distributions from the baseline model. The coupling parameters, fit fractions and other quantities for this model are shown in Table 2. The χ 2 per degree of freedom for this model is calculated to be 40483/32701 = 1.238, which indicates that although this is formally a poor fit, the model is providing a reasonable description of the data given the very large sample size. Three cascade contributions, from a 1 (1260) + , K 1 (1270) − and K(1460) − resonances, are modelled using the three-body running-width treatment described in Sect. 4. The masses and widths of these states are allowed to vary in the fit. The mass, width and coupling parameters for these resonances are presented in Tables 3, 4  Entries/ (0.02 GeV 2 Bands indicate the expectation from the model, with the width of the band indicating the total systematic uncertainty. The total background contribution, which is very low, is shown as a filled area. In figures that involve a single positively-charged pion, one of the two identical pions is selected randomly. Table 2: Fit fractions and coupling parameters for the RS decay D 0 → K − π + π + π − . For each parameter, the first uncertainty is statistical and the second systematic. Couplings g are defined with respect to the coupling to the channel D 0 → [K * (892) 0 ρ(770) 0 ] L=2 . Also given are the χ 2 and the number of degrees of freedom (ν) from the fit and their ratio. on the parametrisation of the running width described by Eq. 5 and of the form factors described by Eq. 6, and thus there is not a straightforward comparison with the values obtained by other experiments. The largest contribution is found to come from the axial vector a 1 (1260) + , which is a result that was also obtained in the Mark III analysis [6]. This decay proceeds via the colour-favoured external W -emission diagram that is expected to dominate this final state.
There are also large contributions from the different orbital angular momentum configurations of the quasi two-body processes D 0 → K * (892) 0 ρ(770) 0 , with a total contribution of around 20%. The polarisation structure of this component is not consistent with naive expectations, with the D wave being the dominant contribution and the overall hierarchy being D > S > P. This result may be compared with that obtained for the study D 0 → ρ(770) 0 ρ(770) 0 in Ref. [31], where the D-wave polarisation of the amplitude was also found to be dominant.
A significant contribution is found from the pseudoscalar state K(1460) − . This resonance is a 2 1 S 0 excitation of the kaon [32]. Evidence for this state has been reported in partial-wave analyses of the process K ± p → K ± π + π − p [33,34], manifesting itself as a 0 − state with mass ≈ 1400 MeV/c 2 and width ≈ 250 MeV/c 2 , coupling to the K * (892) 0 π − and [π − π + ] L=0 K − channels. The intermediate decays of the K(1460) − meson are found to be roughly consistent with previous studies, with approximately equal partial widths to K * (892) 0 π − and [π + π − ] L=0 K − . The resonant nature of this state is confirmed using a model-independent partial-wave analysis (MIPWA), following the method first used by the E791 collaboration [35,36]. The relativistic Breit-Wigner lineshape is replaced by a parametrisation that treats the real and imaginary parts of the amplitude at 15 discrete positions in s K − π + π − as independent pairs of free parameters to be determined by the fit. The amplitude is then modelled elsewhere by interpolating between these values using cubic splines [37]. The Argand diagram for this amplitude is shown in Fig. 3, with points indicating the values determined by the fit, and demonstrates the phase motion expected from a resonance.
Four-body weak decays contain amplitudes that are both even, such as D → [V V ] L=0,2 , where V and V are vector resonances, and odd, such as D → [V V ] L=1 , under parity transformations. Interference between these amplitudes can give rise to parity asymmetries which are different in D 0 and D 0 decays. These asymmetries are the result of strong-phase differences, but can be mistaken for CP asymmetries [38]. Both sources of asymmetry can be studied by examining the distribution of the angle between the decay planes of the two quasi two-body systems, φ, which can be constructed from the three-momenta p of the decay products in the rest frame of the D 0 meson as wheren ab is the direction normal to the decay plane of a two-particle system ab, andp K − π + is the direction of the combined momentum of the K − π + system. The interference between P -even and P -odd amplitudes averages to zero when integrated over the entire phase space. Therefore, the angle φ is studied in regions of phase space. The region of the K * (892) 0 and ρ(770) 0 resonances is studied as the largest P -odd amplitude is the decay D 0 → [K * (892) 0 ρ(770) 0 ] L=1 . Selecting this region allows the identical pions to be distinguished, by one being part of the K * (892) 0 -like system and the other in the ρ(770) 0 -like system. The data in this region are shown in Fig. 4, divided into quadrants of helicity angles, θ A and θ B , defined as the angle between the K − /π − and the D 0 in the rest frame of the K − π + /π − π + system. The distributions show clear asymmetries under reflection about 180 • , indicating parity nonconservation. However, equal and opposite asymmetries are observed in the CP -conjugate mode D 0 → K + π − π − π + , indicating that these asymmetries originate from strong phases, rather than from CP -violating effects. Bands show the expected asymmetries based on the amplitude model, which has been constructed according to the CP -conserving hypothesis, and show reasonable agreement with the data.       0.813 ± 0.032 ± 0.136 112.9 ± 2.6 ± 9.5 β 0 0.315 ± 0.010 ± 0.022 46.7 ± 1.9 ± 3.0  Figure 5: Distributions for six invariant-mass observables in the WS decay D 0 → K + π − π − π + . Bands indicate the expectation from the model, with the width of the band indicating the total systematic uncertainty. The total background contribution is shown as a filled area, with the lower region indicating the expected contribution from mistagged D 0 → K + π − π − π + decays. In figures that involve a single negatively-charged pion, one of the two identical pions is selected randomly.

Results for the WS decay
Invariant mass-squared distributions for D 0 → K + π − π − π + are shown in Fig. 5. Large Table 6: Fit fractions and coupling parameters for the WS decay D 0 → K + π − π − π + . For each parameter, the first uncertainty is statistical and the second systematic. Couplings g are defined with respect to the coupling to the decay D 0 → [K * (892) 0 ρ(770) 0 ] L=2 . Also given are the χ 2 and the number of degrees of freedom (ν) from the fit and their ratio. contributions are clearly seen in s K + π − from the K * (892) 0 resonance. The fit fractions and amplitudes of the baseline model are given in Table 6. The χ 2 per degree of freedom for the fit to the WS data is 350/243 = 1.463. If the true WS amplitude has a comparable structure to the RS amplitude, it contains several decay chains at the O(1%) level that cannot be satisfactorily resolved given the small sample size, and hence the quality of the WS fit is degraded by the absence of these subdominant contributions. Dominant contributions are found from the axial kaons, K 1 (1270) + and K 1 (1400) + , which are related to the same colour-favoured W -emission diagram that dominates the RS decay, where it manifests itself in the a 1 (1260) + K − component. The contribution from the K 1 (1400) + resonance is larger than that from the K 1 (1270) + resonance. It is instructive to consider this behaviour in terms of the quark states, 1 P 1 and 3 P 1 , which mix almost equally to produce the mass eigenstates, where θ K is a mixing angle. The mixing is somewhat less than maximal, with Ref. [39] reporting a preferred value of θ K = (33 +6 −2 ) o . In the WS decay, the axial kaons are produced via a weak current, which is decoupled from the 1 P 1 state in the SU(3) flavour-symmetry limit. If the mixing were maximal, the mass eigenstates would be produced equally, but a smaller mixing angle results in a preference for the K 1 (1400), which is qualitatively consistent with the pattern seen in data. In the RS decay, the axial kaons are not produced by the external weak current, and hence there is no reason to expect either quark state to be preferred. The relatively small contribution from the K 1 (1400) is then understood as a consequence of approximately equal production of the quark states.
The coupling and shape parameters of the K 1 (1270) + resonance are fixed to the values measured in the RS nominal fit. A fit is also performed with these coupling parameters Table 7: Decay chains taken into account in alternative parametrisations of the RS decay mode D 0 → K − π + π + π − . For each chain, the fraction of models in the ensemble that contain this decay, together with the associated average fit fraction, F , are shown. Components are not tabulated if they contribute to all models in the ensemble, or if they contribute to less than 5% of the models.

Decay chain
Fraction F of models [%] [%] free to vary, and the parameters are found to be consistent with those measured in the RS decay. A large contribution is found from D 0 → ρ(1450) 0 K * (892) 0 decays in all models that describe the data well. This contribution resembles a quasi nonresonant component due to the large width of the ρ(1450) 0 resonance, and is likely to be an effective representation of several smaller decay chains involving the K * (892) 0 resonance that cannot be resolved with the current sample size.

Alternative parametrisations
The model-finding procedure outlined in Sect. 5.4 results in ensembles of parametrisations of comparable quality and complexity. The decay chains included in the models discussed above are included in the majority of models of acceptable quality, with further variations made by addition of further small components. The fraction of models in this ensemble containing a given decay mode are shown in Table 7 for the RS decay mode with the average fit fraction associated with each decay chain also tabulated. The ensemble of RS models consists of about 200 models with χ 2 per degree of freedom varying between 1.21 and 1.26. Many of the decay chains in the ensemble include resonances, such as the K 1 (1270) − , decaying via radially excited vector mesons, such as the ρ(1450) 0 and K * (1410) 0 mesons. In particular, the decay K 1 (1270) − → ρ(1450) 0 K − is included in the models discussed in Sects. 6.2 and 6.3 and is found in the majority of the models in the ensemble. This decay channel of the K 1 (1270) − meson has a strong impact at low dipion masses due to the very large width of the ρ(1450) 0 resonance, of about 400 MeV/c 2 . Models excluding this component are presented as alternative parametrisations in Appendix E 22.04 ± 0.28 ± 2.09 21.87 ± 1.51 as this decay mode has not been studied extensively in other production mechanisms of the K 1 (1270) − resonance, and the ensemble contains models without this decay chain of similar fit quality to the baseline model. The situation can be clarified with independent measurements of the properties of these resonances. The a 1 (1640) + resonance is also found in many models in the ensemble, and is likely to be present at some level despite only the low-mass tail of this resonance impacting the phase space. This resonance strongly interferes with the dominant a 1 (1260) + component and, as the parameters of this resonance are poorly known, improved external inputs will be required to correctly constrain this component. The coupling parameters cannot strictly be compared between different models, as in many cases these coupling parameters have a different interpretation depending on the choice of the model. However, it is instructive to consider how the fit fractions vary depending on the choice of model, which is shown in Table 8. It is also useful to consider how the choice of model impacts upon the fitted masses and widths, which is shown in Table 9. The values for the model described in Sect. 6.2 are also shown, which has compatible values with the ensemble. The variation with respect to the choice of model is characterised by the RMS of the parameters in the ensemble, and is of a comparable size to the combined systematic uncertainty from other sources on these parameters.
The D 0 → K + π − π − π + ensemble consists of 108 models, all of which have a χ 2 per degree of freedom of less than 1.45, with the best models in the ensemble having a χ 2 per degree of freedom of about 1.35. The fraction of models in this ensemble containing a given decay mode are shown in Table 10. In particular, there should be percent-level contributions from some of the decay chains present in the D 0 → K − π + π + π − mode, such as D 0 → a 1 (1260) − K + and D 0 → K * (892) 0 [π + π − ] L=0 . In addition to the marginal decays of the K 1 (1270) + present in the D 0 → K + π − π − π + ensemble, the models suggest contributions from the K * (1680), which resembles a nonresonant component due to its large width and position on the edge of the phase space. As is the case for the large D 0 → K * (892) 0 ρ(1450) component, this contribution is likely to be mimicking several smaller decay channels that cannot be resolved with the current sample size. Table 10: Decay chains taken into account in alternative parametrisations of the WS decay mode D 0 → K + π − π − π + . For each chain, the fraction of models in the ensemble that contain this decay, together with the associated average fit fraction, F , are shown. Components are not tabulated if they contribute to all models in the ensemble, or if they contribute to less than 5% of the models.

Coherence factor
The coherence factor R K3π and average strong-phase difference δ K3π are measures of the phase-space-averaged interference properties between suppressed and favoured amplitudes, and are defined as [40] R where A(D 0 → K + 3π) is the amplitude of the suppressed decay and A(D 0 → K + 3π) is the favoured amplitude for D 0 decays. Additionally, it is useful to define the average ratio of amplitudes as Knowledge of these parameters is necessary when making use of this decay in B − → DK − transitions for measuring the CP -violating phase γ [40], and can also be exploited for charm mixing studies. Observables with direct sensitivity to the coherence factor and related parameters have been measured in e + e − collisions at the ψ(3770) resonance with CLEO-c data [41], and through charm mixing at LHCb [4]. A global analysis of these results [41] yields The baseline models presented in Sect. 6 can be used to calculate the model-derived coherence factor R mod K3π = 0.458 ± 0.010 ± 0.012 ± 0.020, where the first uncertainty is statistical, the second systematic, and the third the uncertainty from the choice of WS model. This uncertainty is assigned by taking the spread in values from an ensemble of alternative models from the model-building algorithm, requiring that models have a χ 2 per degree of freedom of less than 1.5, and that all unconstrained components in the fit have a significance of > 2σ. This result is in good agreement with the direct measurement in Ref. [41]. This analysis has no sensitivity to δ K3π and r K3π as each amplitude contains an arbitrary independent amplitude and phase.
The stability of the local phase description can also be verified by evaluating the modelderived coherence factor and associated parameters in different regions of phase space. This is equivalent to changing the definition of Eq. 19 such that integrals are performed over a limited region of phase space. In this case, it is also possible to determine the local values of δ K3π and r K3π relative to the phase-space averaged values. Therefore, overall normalisation factors are fixed such that the central values of the direct measurement are correctly reproduced.
In order to define these regions, the phase space is divided into hypercubes using the algorithm described in Sect. 5.2. The division is done such that the hypercubes cannot be smaller in any dimension than 50 MeV/c 2 . The hypercubes are grouped into bins of average phase difference between the two amplitudes in the bin, using the amplitude models described in Sect. 6. The range [−180 • , 180 • ] in phase difference between the two decay modes is split into eight bins. The division of this range is done such that each bin is expected to have an approximately equal population of WS events within the bin. The coherence factors, average strong phases and amplitude ratios and their RMS spread arising from the choice of WS model are summarised in Table 11. Good stability with respect to the choice of model is observed, which is a consequence of the dominant features of the amplitude being common for all models, and gives confidence to using the models presented in this paper to define regions of interest for future binned measurements of γ or studies of charm mixing. The relatively high coherence factor in some regions of phase-space demonstrates the potential improvements in sensitivity to measurements of CP -violating observables.

Conclusions
The four-body decay modes D 0 → K ∓ π ± π ± π ∓ have been studied using high-purity timeintegrated samples obtained from doubly tagged B → D * + (2010)[D 0 π + ]µX decays. For the RS decay mode D 0 → K − π + π + π − , the analysis is performed with a sample around sixty times larger than that exploited in any previous analysis of this decay. For the WS mode D 0 → K + π − π − π + , the resonance substructure is studied for the first time with ≈ 3000 signal candidates. Both amplitude models are found to have large contributions from axial resonances, the decays D 0 → a 1 (1260) + K − and D 0 → K 1 (1270/1400) + π − for D 0 → K − π + π + π − and D 0 → K + π − π − π + , respectively. This is consistent with the general picture that W -emission topologies are crucial in describing these decays. Interference between the parity-even and parity-odd amplitudes causes large local parity violations, which are shown to be reasonably well modelled in the RS decay. A significant contribution from the pseudoscalar resonance K(1460) − is identified, which is validated using the modelindependent partial waves method.
The coherence factor is calculated using the models, and is found to be consistent with direct measurements. It is found that the calculated value is relatively stable with respect to the parametrisation of subdominant amplitudes in the WS model. These models therefore provide a valuable input to future binned measurements of the CP -violating parameter γ and charm-mixing studies.

Appendices A Spin formalism
The effects of spin and orbital angular momentum are calculated using the Rarita-Schwinger formalism, following a similar prescription to that described in Ref. [23]. Spin-matrix elements for quasi two-body processes are constructed in terms of a series of polarisation and pure orbital angular momentum tensors. Consider the decay of particle a that has integer spin J, into particles b and c, which have integer spin s b , s c , respectively. All three particles have an associated polarisation tensor, (a,b,c) , which is of rank equal to the spin of the particle. The decay products b, c will also in general have a relative orbital angular momentum l, which is expressed in terms of the pure orbital angular momentum tensor, L µ...ν , which is of rank l. The matrix element for the decay is where the tensor G ... combines the polarisation and pure orbital angular momentum tensor to produce a scalar object. This tensor is constructed from combinations of the metric tensor g µν and the Levi-Civita tensor contracted with the four-momenta of the decaying particle, ε µναβ P µ . The second of these tensors is used only if J − (l − s b − s c ) is odd, and ensures that matrix elements have the correct properties under parity transformations. The matrix element can also be written by defining the current, J , of the decaying particle: where the µ represents a set of Lorentz indices µ...ν, a shorthand which will be used throughout this section. The isobar model factorises an N -body decay into a sequence of two-body processes. Each of these quasi two-body decays can be described with a single spin matrix element, and hence the total matrix element is the product of N − 1 matrix elements: For example, consider the quasi two-body decay P → X[ab]Y [cd]. The matrix element for this decay is where the sums are over the possible polarisations of the intermediate states.
It is preferable to build a generic formulation of the total matrix element for arbitrary topologies, spins and angular momenta, rather than performing an explicit computation for each possible process. A generic approach to computing matrix elements is to introduce a generalised "current" associated with a decaying particle that has absorbed the matrix elements of its decay products. This current can be written in terms of the currents of its decay products as where S 1,2 µ is the spin-projection operator of decay products (1,2), which has been used to sum intermediate polarisation tensors, using the definition The first few projection operations, which are sufficient for describing charm decays, are This operator projects out the component of a tensor that is orthogonal to the fourmomentum of a particle, and has rank 2J for a particle of spin J. The orbital angular momentum tensors are also constructed from the spin projection operators and the relative momentum of the decay products, Q a [23], and are written as The matrix element for a generic cascade of particle decays can then be calculated recursively. In the case of the decay of a spinless particle, the matrix element for the total decay process is identical to the current of the decaying particle. The generalised current is therefore merely a convenient device for organising the computation of spin matrix elements, but is not in general associated with the propagation of angular momentum. It is also useful to define the spin-projected currents, S µν J ν , which will be written as S, V µ , T µν for (pseudo)scalar, (pseudo)vector and (pseudo)tensor states, respectively.
The rules for how the different spin-projected currents are written in terms of each other is given in Table 12, where these relations are derived by considering the symmetries of Lorentz indices and the parity properties of the matrix element. All of the coupling structures necessary to describe P → 4P are uniquely determined by these constraints, although this property does not hold in general. This allows complicated spin configurations to be calculated in terms of a simple and consistent set of rules. The rules are written with consistent dependencies to clarify their derivations, and in some cases simplified forms are also given. These simplifications typically rely on the symmetry properties of the Levi-Civita tensor and the relationship S ab S bc = S a c , which is the defining characteristic of a projection operator. Table 12: Rules for calculating the current associated with a given decay chain in terms of the currents of the decay products. Where relevant, the spin projection operator S and the orbital angular momentum operators L are those for the decaying particle.

Topology
Current Simplified current

B List of decay chains
The list of possible decay chains is built from what is allowed by the relevant conservation laws. Approximately one hundred different decay chains modes are included as possible contributions to the model. Certain cascade decays already have well known sub-branching ratios. For example, although the K 1 (1400) decays almost exclusively via the K * (892), the various decays of the K 1 (1400) are treated separately without assumption about their branching ratios. The different components can be split into the same groups as in Sect. 4: where Y ππ is one of the following states: ρ(770) 0 , ρ(1450) 0 , f 2 (1270) or [π + π − ] L=0 , and Y Kπ is one of the following: K * (892) 0 , K * (1410) 0 , The [π + π − ] L=0 and [K ∓ π ± ] L=0 contributions are modelled using K matrices. In cases with a scalar contribution and a radial recurrence of a vector state, such as ρ(1450) 0 [K ∓ π ± ] L=0 , the K matrix is fixed to be the same as the first vector, i.e. the K-matrix parameters of ρ(770) 0 [K ∓ π ± ] L=0 . For vector-vector and vector-tensor contributions, the different possible polarisation states are included together in the model building. The contributions from the radial excitations of the kaon are only included as a possibility when included with the π + π − S-wave, as the other decay chains involving this resonance, for example the decay K * (1410) 0 ρ(770) 0 , tend to have large interference terms, which requires fine tuning with other amplitudes and hence are considered to be unphysical.
where X Kππ is one of the following states: All of these states are considered under all possible orbital configurations that obey the respective conservation laws. The various contributions assigned for different systematic uncertainties are summarised in this appendix by a series of tables. The legend for these is given in Table 13, including which sources of uncertainty are considered on each decay mode. The breakdown of systematic uncertainties for the RS decay D 0 → K − π + π + π − for coupling parameters, fit fractions and other parameters are given in Tables 14 and 15 for the quasi two-body decay chains and cascade decay chains, respectively. The systematic uncertainties for the WS mode D 0 → K + π − π − π + are given in Table 16 for both coupling parameters and the fit fractions.

D Interference fractions
The interference fraction between decay chains a and b is where the sum over j is over all of the decay chains. For cascade processes, the different secondary isobars contribute coherently to the interference fractions. The interference fractions are presented in Tables 17 and 18 for RS and WS decay modes, respectively. For each decay mode, the largest interference fractions are between the axial vector decay chain, and the lowest orbital angular momentum vector-vector decay chain.