Amplitude analyses of D0 → π+π−π+π− and D0 → K+K−π+π− decays

The resonant substructure of D0 → π+π−π+π− decays is studied using data collected by the CLEO-c detector. An amplitude analysis is performed in order to disentangle the various intermediate state contributions. To limit the model complexity a data driven regularization procedure is applied. The prominent contributions are the decay modes D0 → a1(1260)+ π−, D0 → σ f0(1370) and D0 → ρ(770)0ρ(770)0. The broad resonances a1(1260)+, π(1300)+ and a1(1640)+ are studied in detail, including quasi-modelindependent parametrizations of their lineshapes. The mass and width of the a1(1260)+ meson are determined to be ma1(1260)+ = [1225 ± 9 (stat) ± 17 (syst) ± 10 (model)] MeV/c2 and Γa1(1260)+ = [430 ± 24 (stat) ± 25 (syst) ± 18 (model)] MeV. The amplitude model of D0 → K+K−π+π− decays obtained from CLEO II.V, CLEO III, and CLEO-c data is revisited with improved lineshape parametrizations. The largest components are the decay modes D0 → ϕ(1020)ρ(770)0, D0 → K1(1270)+K− and D0 → K(1400)+K−. The fractional CP -even content of the decay D0 → π+π−π+π− is calculated from the amplitude model to be F+4π = [72.9 ± 0.9(stat) ± 1.5(syst) ± 1.0(model)] %, consistent with that obtained from a previous model-independent measurement. For D0 → K+K−π+π− decays, the CP -even fraction is measured for the first time and found to be F+KKππ = [75.3 ± 1.8 (stat) ± 3.3 (syst) ± 3.5 (model)] %. The global decay rate asymmetries between D0 and D¯0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\overline{D}}^0 $$\end{document} decays are measured to be ACP4π=+0.54±1.04stat±0.51syst%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathcal{A}}_{CP}^{4\uppi}=\left[+0.54\pm 1.04\ \left(\mathrm{stat}\right)\pm 0.51\ \left(\mathrm{syst}\right)\right]\% $$\end{document} and ACPKKππ=+1.84±1.74stat±0.30syst%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathcal{A}}_{CP}^{KK\pi \pi}=\left[+1.84\pm 1.74\ \left(\mathrm{stat}\right)\pm 0.30\ \left(\mathrm{syst}\right)\right]\% $$\end{document}. A search for CP asymmetries in the amplitude components yields no evidence for CP violation in either decay mode.


Introduction
We present amplitude analyses for D 0 → h + h − π + π − decays, where h ± is either a pion or a kaon. These decay modes have the potential to make an important contribution to the determination of the CP -violating phase γ (φ 3 ) ≡ − arg(V ud V * ub /V cd V * cb ) in B − → DK − and related decays [1][2][3][4][5][6]. The all-charged final states (impossible in three-body decays of D 0 ) particularly suit the environment of hadron collider experiments, such as LHCb. The sensitivity to the weak phase can be significantly improved with a measured D-decay amplitude model, either to be used directly in the γ extraction, or in order to optimize model-independent measurements [4,[7][8][9][10].
A study of the rich resonance structure of these four-body decays is also of considerable interest in its own right. Figure 1 shows the dominant processes that contribute to the visible structure in the phase space. The color-favored tree diagram manifests as a cascade whereby a resonance decays into another resonance before decaying into the final state. Due to the identical quark content produced in the weak and spectator interactions, a given process and its CP -conjugate may arise even from the same initial state. Such processes, which we refer to as non-self-conjugate, are also known as flavor-non-specific decays as flavor-tagging is required to distinguish between the source of these two partners despite not being CP eigenstates. The color-suppressed tree diagram and the W -exchange diagram result in self-conjugate intermediate states such as ρ(770) 0 ρ(770) 0 or ρ(770) 0 φ(1020) whose partial waves are eigenstates of CP . Certain intermediate states in D 0 → K + K − π + π − decays, for instance K * (892) 0K * (892) 0 , are only accessible via the W -exchange diagram.
The decay D 0 → π + π − π + π − provides an excellent environment to study the properties of the a 1 (1260) + meson, whose width is an unresolved question, currently given as 250 − 600 MeV in the Particle Data Group's Review of Particle Physics (PDG) [11]. The only previous analysis of the D 0 → π + π − π + π − amplitude structure was published by the FOCUS collaboration based on approximately 6000 D 0 , D 0 → π + π − π + π − signal events [12]. The analysis presented here benefits from the ability to distinguish D 0 from D 0 decays and a larger data sample of approximately 7000 signal events.
Based on the four-body amplitude formalism and analysis software used in the D 0 → K + K − π + π − amplitude analysis performed by the CLEO collaboration [13], we introduce significant improvements especially in the parametrization of three-body resonances. Using a state-of-the-art parametrization of the a 1 (1260) + lineshape, we present new measurements of the a 1 (1260) + mass and width. By utilizing different parametrizations, we confirm a significant dependence of the measured width on the lineshape itself. We also observe contributions from the decay modes D 0 → a 1 (1640) + π − and D 0 → π(1300) + π − , not seen in previous analyses and provide model-independent complex lineshapes for the a 1 (1260) + , a 1 (1640) + and π(1300) + mesons.
In addition to our new D 0 → π + π − π + π − analysis, we also revisit the CLEO D 0 → K + K − π + π − data using the improved formalism and analysis procedures presented in this paper. Prior to the CLEO analysis, an amplitude analysis of the decay D 0 → K + K − π + π − was also performed by FOCUS [14].  Figure 1. Examples of the color-favored (a), color-suppressed (b) and W -exchange (c) diagrams that contribute towards the resonant structure in D 0 → π + π − π + π − and D 0 → K + K − π + π − decays.
This article is structured as follows: after an introduction to the CLEO II.V, CLEO III, and CLEO-c experiments in section 2 and a description of the event selection in section 3, the amplitude formalism and its implementation is described in section 4 and section 5. The results of the fit to data, including a model-dependent measurement of the fractional CP -even content and search for direct CP violation, are presented in section 6 and section 7. Systematic uncertainties are outlined in section 8, and our conclusions are given in section 9. Additional technical details of the analyses can be found in the appendices and supplementary material.

Data set and CLEO detector
The data analyzed in this paper were produced in symmetric e + e − collisions at CESR between 1995 and 2008, and collected with three different configurations of the CLEO detector: CLEO II.V, CLEO III, and CLEO-c.
In CLEO II.V [15,16] tracking was provided by a three-layer double-sided silicon vertex detector, and two drift chambers. Charged particle identification came from dE/dx information in the drift chambers, and time-of-flight (TOF) counters inserted before the calorimeter. For CLEO III [17] a new silicon vertex detector was installed, and a ring imaging Cherenkov (RICH) detector was deployed to enhance the particle identification abilities [18]. In CLEO-c, the vertex detector was replaced with a low-mass wire drift chamber [19]. A superconducting solenoid supplied a 1.5 T magnetic field for CLEO II.V and III, and 1 T for CLEO-c operation, where the average particle momentum was lower. In all detector configurations, neutral pion and photon identification was provided by a 7800-crystal CsI electromagnetic calorimeter.
Four distinct data sets are analyzed in the present study:

JHEP05(2017)143
(1) approximately 9 fb −1 accumulated at √ s ≈ 10 GeV by the CLEO II.V detector; (2) a total of 15.3 fb −1 accumulated by the CLEO III detector in an energy range √ s = 7.0 − 11.2 GeV, with over 90% of this sample taken at √ s = 9.5 − 10.6 GeV; (3) 818 pb −1 collected at the ψ(3770) resonance by the CLEO-c detector; (4) a further 600 pb −1 taken by CLEO-c at √ s = 4170 MeV, where √ s is the total energy delivered by the beam in the center-of-mass system (CMS). These samples are referred to as the CLEO II.V, CLEO III, CLEO-c 3770 and CLEO-c 4170 data sets, respectively.
Detector response is studied with GEANT-based [20] Monte Carlo (MC) simulations of each detector configuration, in which the MC events are processed with the same reconstruction algorithm as used for data.

Event selection
We select events where one neutral D meson decays either into a π + π − π + π − or K + K − π + π − final state. The analysis considers two classes of signal decays, for both of which information on the quantum numbers of the meson decaying to the signal mode is provided by an event tag.
(i) Flavor-tagged decays are selected from the CLEO II.V and CLEO III data sets, in which the flavor of the decaying meson is determined by the charge of the 'slow pion', π s , in the D * + → D 0 π + s decay chain. Flavor-tagged decays are also selected from the two CLEO-c data sets, where here the tag is obtained through the charge of a kaon associated with the decay of the other D meson in the event. The wrong tag fractions for each data set are represented by the parameter w, given in ref. [13].
(ii) CP -tagged decays are selected in the CLEO-c 3770 data set alone. In ψ(3770) decays the D − D pair is produced coherently. Therefore, the CP of the signal D can be determined if the other D meson is reconstructed in a decay to a CP -eigenstate. Useful information is also obtained if the tagging meson is reconstructed decaying into the modes K 0 S π + π − or K 0 L π + π − , for which the relative contribution of CP -even and CP -odd states is known [21].
The D 0 → π + π − π + π − analysis uses only the flavor-tagged subset of the CLEO-c 3770 data sample, while D 0 → K + K − π + π − makes use of all the data sets described. The selection criteria for producing the data sets of each of these classes is discussed in detail in ref. [13] and is identical to that used in our analysis, except for a few improvements that will be highlighted where applicable.

JHEP05(2017)143
incoherent process since the K 0 S lifetime is much longer than those of any other possible intermediate resonance. Therefore, K 0 S decays are rejected if the invariant mass of any π + π − combination is within 7.5 MeV/c 2 of the world-average K 0 S mass [11]. Two nearly uncorrelated kinematic variables are used to define a signal and two sideband background regions. These variables are defined as the beam-constrained mass, where p D is the reconstructed three-momentum of the candidate D in the CMS; and the missing energy ∆E, where E D is the total reconstructed energy of candidate D in the CMS. Signal events should have missing energy close to zero and beam-constrained mass close to that of the nominal D 0 mass, m D [11]. By construction, the m bc width is a measure of the beamenergy spread while the ∆E width is dominated by the detector resolution. Candidates that satisfy m bc > 1.83 GeV/c 2 and |∆E| < 0.1 GeV are retained for further analysis.
As the sideband events are used to study the background contribution within the signal region, it is crucial to select signal and background regions with a mutual and constant invariant mass, i.e. that of the D meson. First, a region of constant invariant mass is obtained by selecting events with This relation describes an annulus in m bc and ∆E space. Lines normal to this annulus of constant invariant mass have an angle of inclination about the center of the annulus. A signal region around the D mass peak is then defined by requiring |θ − θ D | < 0.004, where m bc = m D and ∆E = 0 GeV at θ D , as shown in figure 2.
Similarly, sideband regions are defined with |θ − θ D | > 0.006. These criteria preserve the range of invariant mass selected throughout the kinematic variables m bc and ∆E, ensuring the distribution of events in phase space are consistent between regions. The signal region contains 9247 D → π + π − π + π − candidates. To estimate the signal purity of the sample, a two-dimensional unbinned maximum likelihood fit to m bc and ∆E is performed in the whole range. While the signal peak is modeled with a sum of three (two) Gaussian functions, the combinatorial background is described by an ARGUS [22] (linear) function in m bc (∆E). The number of signal events within the signal region is estimated from the fit result displayed in figure 3, to be 7250 ± 56 (stat) ± 46 (syst) events, where the first uncertainty is statistical and second is systematic. The signal fraction f Sig , in this region is f Sig = (78.4±0.6 (stat)±0.5 (syst))%. These systematic uncertainties are estimated by repeating the fit with different appropriate probability density function (PDF) hypotheses. As we observed a negligible impact of the background on our analysis, further improvements of the signal purity were not studied. . Distribution of D 0 → π + π − π + π − candidate events in missing energy ∆E and beam constrained mass within the selection regions, which are bounded by the annulus of constant invariant mass and lines normal to it. The central region (blue) is defined as the signal region, with sideband regions (red) providing background samples.  . Beam constrained mass (a) and missing energy (b) distribution of D 0 → π + π − π + π − candidates, overlaid with the projections of the fitted PDF (solid black line). The signal component is shown in blue (dashed) and the background component in red (dashed).

D
With respect to ref. [13], we veto the π + π − invariant mass region around the K 0 S mass, which removes essentially all peaking background from D 0 → K 0 S (→ π + π − ) K + K − , greatly simplifying our analysis. The K 0 S veto depends on the CLEO configuration, as the mass resolution is better for data collected with the CLEO-c configurations. For data collected with CLEO (CLEO-c), the π + π − invariant mass combination does not fall within 16.5 (12) MeV/c 2 of the world-average K 0 S mass. In addition, for the flavor-tagged data, several changes have been applied with respect to ref. [13]. The CLEO II.V minimum track momenta cut for the D daughters is raised Signal region Sideband CLEO II.V 144.6 MeV/c 2 < ∆m < 146.2 MeV/c 2 148.5 MeV/c 2 < ∆m < 160.0 MeV/c 2 CLEO III 144.6 MeV/c 2 < ∆m < 146.1 MeV/c 2 148.5 MeV/c 2 < ∆m < 160.0 MeV/c 2 CLEO-c 3770 |m bc − m D | < 0.005 GeV/c 2 1.834 GeV/c 2 < m bc < 1.854 GeV/c 2 1.876 GeV/c 2 < m bc < 1.890 GeV/c 2 CLEO-c 4170 2.005 GeV/c 2 < m bc < 2.030 GeV/c 2 1.880 GeV/c 2 < m bc < 1.920 GeV/c 2 Table 1. Signal region and sideband definitions in the ∆m or m bc kinematic variable, for flavortagged D 0 → K + K − π + π − data in the different CLEO configurations.

Sample
Signal  Table 2. Updated number of signal candidates and fractions in the signal region, for flavor-tagged to 275 MeV/c as the MC was found not to represent the data sufficiently well below this value. As in ref. [13], the kinematic variables that describe signal in the CLEO II.V and CLEO III samples are the reconstructed D mass m KKππ , and the mass difference between the D * and D candidates, ∆m. We take advantage of the possibility to ensure a constant  [23]. The procedure to measure the purity in the signal region of each sample is identical to that of the previous analysis [13]. The events retained for the amplitude analysis and signal fractions for the improved selection criteria are given in table 2.

JHEP05(2017)143
Key differences are in the formalism used for the spin factors, where we now use a more consistent and intuitive implementation of the Zemach formalism [26][27][28], and an improved description of the lineshapes of resonances decaying to three-body final states.
The differential decay rate of a D 0 meson with mass, m D 0 , decaying into four pseudoscalar particles with four-momenta p i = (E i , p i ) (i = 1, 2, 3, 4) is given by where the transition amplitude A D 0 (x), describes the dynamics of the interaction, dΦ 4 is the four-body phase space element [29], and x represents a unique set of kinematic conditions within the phase space of the decay. Each final state particle contributes three observables, manifesting in their three-momentum, summing up to twelve observables in total. Four of them are redundant due to four-momentum conservation and the overall orientation of the system can be integrated out. The remaining five independent degrees of freedom unambiguously determine the kinematics of the decay. Convenient choices for the kinematic observables include the invariant mass combinations of the final state particles or acoplanarity and helicity angles [30,31]. It is however important to take into account that, while m 2 12 , m 2 23 are sufficient to fully describe a three-body decay, the obvious extension to four-body decays with m 2 ij , m 2 ijk requires additional care, as these variables alone are insufficient to describe the parity-odd moments possible in four-body kinematics.
In practice, we do not need to choose a particular five-dimensional basis, but use the full four-vectors of the decay in our analysis. The dimensionality is handled by the phase space element which can be written in terms of any set of five independent kinematic observables, x = (x 1 , . . . , x 5 ), as is the phase space density. In contrast to three-body decays, the four-body phase space density function is not flat in the usual kinematic variables. Therefore, an analytic expression for φ 4 is taken from ref. [32].
The total amplitude for the D 0 → h 1 h 2 h 3 h 4 decay is given by the coherent sum over all intermediate state amplitudes A i (x), each weighted by a complex coefficient a i = |a i | e i φ i to be measured from data, To construct A i (x), the isobar approach is used, which assumes that the decay process can be factorized into subsequent two-body decay amplitudes [33][34][35]. This gives rise to two different decay topologies; quasi two-body decays In either case, the intermediate state amplitude is parameterized as a product of form factors B L , included for each vertex of the -8 -

JHEP05(2017)143
decay tree, Breit-Wigner propagators T R , included for each resonance R, and an overall angular distribution represented by a spin factor S, (4.5) As the π + π − π − π + final state involves two pairs of indistinguishable pions, the amplitudes are Bose-symmetrized and therefore symmetric under exchange of like-sign pions. We define the CP -conjugate phase space point x such that it is mapped onto x by the interchange of final state charges, and the reversal of three-momenta. If x, x are expressed as a function of the four-momenta (E i , p i ) (where i labels the particle), this implies for and equivalently for D 0 → π + π − π + π − . The CP -conjugate of a given intermediate state and the total D 0 decay amplitude is defined as Unless stated otherwise, we assume CP conservation in the D 0 decay, implyingā i = a i . Moreover, CP conservation in the strong interaction is implemented in the cascade topology by the sharing of couplings between related quasi-two-body final states. For example, given the two a i parameters required for D 0 → π − a 1 (1260) + with a 1 (1260) + → ρ(770) 0 π + and a 1 (1260) + → σ π + , the amplitude D 0 → π + a 1 (1260) − with a 1 (1260) − → ρ(770) 0 π − and a 1 (1260) − → σ π − only requires one additional global complex parameter to represent the different weak processes of D 0 → a 1 (1260) + π − and D 0 → a 1 (1260) − π + , while the relative magnitude and phase of a 1 (1260) − → ρ(770) 0 π − and a 1 (1260) − → σ π − are the same as for a 1 (1260) + → ρ(770) 0 π + and a 1 (1260) + → σ π + . For historical reasons, this constraint is only applied to the π + π − π + π − final state, but, as discussed in section 7, the results we obtain for the K + K − π + π − final state are also compatible with CP conservation in the strong interaction.

Form factors and resonance lineshapes
To account for the finite size of the decaying resonances, the Blatt-Weisskopf penetration factors, derived in ref. [36] by assuming a square well interaction potential with radius r BW , are used as form factors, B L . They depend on the breakup momentum q, and the orbital angular momentum L, between the resonance daughters. Their explicit expressions are B 0 (q) = 1, B 2 (q) = 1/ 9 + 3 (q r BW ) 2 + (q r BW ) 4 . Resonance lineshapes are described as function of the energy-squared, s, by Breit-Wigner propagators , (4.10) featuring the energy-dependent mass M (s) (defined below), and total width, Γ(s). The latter is normalized to give the nominal width, Γ 0 , when evaluated at the nominal mass m 0 , i.e. Γ 0 = Γ(s = m 2 0 ). For a decay into two stable particles R → AB, the energy dependence of the decay width can be described by where q 0 is the value of the breakup momentum at the resonance pole [37]. The energy-dependent width for a three-body decay R → ABC, on the other hand, is considerably more complicated and has no analytic expression in general. However, it can be obtained numerically by integrating the transition amplitude-squared over the phase space, and therefore requires knowledge of the resonant substructure. The three-body amplitude A R→ABC can be parameterized similarly to the four-body amplitude in eq. (4.5). In particular, it includes form factors and propagators of intermediate two-body resonances. Both eq. (4.11) and eq. (4.12) give only the partial width for the decay into a specific channel. To obtain the total width, a sum over all possible decay channels has to be performed, Γ(s) = i g i Γ i (s), (4.13) where the coupling strength to channel i, is given by g i . Branching fractions B i are related to the couplings g i via the equation [11] (4.14) As experimental values are usually only available for the branching fractions, eq. (4.14) needs to be inverted to obtain values for the couplings. In practice, this is solved by minimizing the quantity where I i (g) denotes the righthand side of eq. (4.14).
The energy-dependent mass follows from the decay width via the Kramers-Kronig dispersion relation [38,39]: Here, the energy-dependent mass is normalized such that M 2 (s = m 2 0 ) = m 2 0 . In practice, the energy-dependent mass is often approximated as being constant, i.e. M 2 (s) = m 2 0 , since -10 -

JHEP05(2017)143
its calculation requires a detailed understanding of the decay width for arbitrarily large energies and is computationally expensive. This is usually justified as the energy-dependent mass needs to satisfy the condition, dM 2 (s) ds such that M 2 (s) is indeed, approximately constant near the on-shell mass [40]. Larger dispersive effects are thus only expected for very broad resonances. The treatment of the lineshape for various resonances considered in this analysis is described in what follows. The nominal masses and widths of the resonances are taken from the PDG [11] with the exceptions described below. We assume an energy-independent mass unless otherwise stated.
For the broad scalar resonance σ, the model from Bugg is used [41]. Besides σ → ππ decays, it includes contributions from the decay modes σ → KK, σ → ηη and σ → ππππ as well as dispersive effects due to the channel opening of the latter. We use the Gournaris-Sakurai parametrization for the ρ(770) 0 → ππ propagator which provides an analytical description of the dispersive term, M 2 (s) [42]. The energy-dependent width of the f 0 (980) resonance is given by the sum of the partial widths into the ππ and KK channels [43], (4.17) where the coupling constants g ππ and g KK , as well as the mass and width are taken from a measurement performed by the BES Collaboration [44]. The total decay widths for both the f 2 (1270) and the f 0 (1370) meson take the channels ππ, KK, ηη and ππππ into account. While the two-body partial widths are described by eq. (4.11), a model for the partial width for a decay into four pions is taken from ref. [45]. The corresponding branching fractions are taken from the PDG [11]. The nominal mass and width of the f 0 (1370) resonance are taken from an LHCb measurement [46]. Equation (4.11) is used for all other resonances decaying into a two-body final state.
To describe the decay width of the axial vector resonance a 1 (1260), the decay channels πππ and KKπ are considered, a 1 (1260)→KKπ (s), (4.18) where isospin symmetry is assumed, i.e. Γ a 1 (1260) + →π 0 π 0 π + (s). The partial width Γ (3) a 1 (1260)→KKπ (s) is calculated from eq. (4.12) assuming the decay proceeds entirely via a 1 (1260) → K * (892) K. The corresponding branching fraction is taken from a CLEO analysis of hadronic τ decays [47]. The calculation of the partial width Γ (3) a 1 (1260)→πππ (s) is more complicated due to the fact that it requires information about the three pion Dalitz plot structure of the a 1 (1260) resonance whose determination in turn, needs the propagator as input. For this reason, we follow an iterative approach. The initial amplitude fit, described in section 6, is performed using an energy-dependent width distribution derived from an uniform phase space population. Afterwards, the energy-dependent  width is recalculated with the results of the substructure analysis and the amplitude fit is subsequently repeated with the new propagator. It is found that the energy-dependent width is not highly sensitive to the details of the Dalitz plot as this procedure converges after a few iterations. As the a 1 (1260) resonance is very broad, the dispersive term is calculated as well. Figure 4 shows the final iteration of the energy-dependent width and mass. The energy-dependent width varies strongly around s ≈ 0.8 GeV 2 where the energy of the π + π − subsystem is equal to the ρ(770) 0 on-shell mass. Around s = 2 GeV 2 , a small hump develops due to the opening of the KKπ channel. The energy-dependent mass indeed shows a plateau around the nominal mass as expected. Note that as the condition of eq. (4.16) is not explicitly enforced by eq. (4.15), it serves as an independent check of whether the main thresholds have been included [38,47].
For the resonances π(1300), a 1 (1640) and π 2 (1670), the energy-dependent width is obtained via the same iterative procedure as for the a 1 (1260) resonance. In case of the π 2 (1670) meson, the KKπ and ωρ(770) 0 thresholds are included with the PDG branching fractions taken from ref. [11], otherwise only decays to three pions are considered. In the D 0 → K + K − π + π − analysis, resonant decays of the K 1 (1270) and K 1 (1400) mesons into the Kρ(770) 0 , K * (892)π, K * 0 (1430)π, Kf 0 (1370) and Kω decay channels are taken into account assuming the lowest possible angular momentum state. For the purpose of evaluating the energy-dependent widths of the excited kaons, these decay channels are assumed to be incoherent and the branching fractions from the PDG are used [11]. The same procedure is applied to obtain the energy-dependent width for the K * (1410) and K * (1680) resonances. In their case, the decay channels Kρ(770) 0 , K * (892)π and Kπ are considered. Some particles may not originate from a resonance but are in a state of relative orbital angular momentum. We denote such non-resonant states by surrounding the particle -12 -JHEP05(2017)143 system with brackets and indicate the partial wave state with an subscript; for example (ππ) S refers to a non-resonant di-pion S-wave. The lineshape for non-resonant states is set to unity.

Spin densities
The spin amplitudes are phenomenological descriptions of decay processes that are required to be Lorentz invariant, compatible with angular momentum conservation and, where appropriate, parity conservation. They are constructed in the covariant Zemach (Rarita-Schwinger) tensor formalism [26][27][28]. At this point, we briefly introduce the fundamental objects of the covariant tensor formalism which connect the particle's four-momenta to the spin dynamics of the reaction and give a general recipe to calculate the spin factors for arbitrary decay trees. Further details can be found in refs. [48,49].
A spin-S particle with four-momentum p, and spin projection λ, is represented by the polarization tensor (S) (p, λ), which is symmetric, traceless and orthogonal to p. These so-called Rarita-Schwinger conditions reduce the a priori 4 S elements of the rank-S tensor to 2S + 1 independent elements in accordance with the number of degrees of freedom of a spin-S state [27,50].
The spin projection operator , for a resonance R, with spin S = {0, 1, 2}, and four-momentum p R , is given by [49]: where g µν is the Minkowski metric. Contracted with an arbitrary tensor, the projection operator selects the part of the tensor which satisfies the Rarita-Schwinger conditions. For a decay process R → AB, with relative orbital angular momentum L, between particle A and B, the angular momentum tensor is obtained by projecting the rank-L tensor q ν 1 R q ν 2 R . . . q ν L R , constructed from the relative momenta q R = p A − p B , onto the spin-L subspace, Their | q R | L dependence accounts for the influence of the centrifugal barrier on the transition amplitudes. For the sake of brevity, the following notation is introduced, Following the isobar approach, a four-body decay amplitude is described as a product of two-body decay amplitudes. Each sequential two-body decay R → A B, with relative -13 -

JHEP05(2017)143
orbital angular momentum L AB , and total intrinsic spin S AB , contributes a term to the overall spin factor given by Here, a polarization vector is assigned to the decaying particle and the complex conjugate vectors for each decay product. The spin and orbital angular momentum couplings are described by the tensors P (S AB ) (R) and L (L AB ) (R), respectively. Firstly, the two spins S A and S B , are coupled to a total spin-S AB state, Φ(x|S AB ), by projecting the corresponding polarization vectors onto the spin-S AB subspace transverse to the momentum of the decaying particle. Afterwards, the spin and orbital angular momentum tensors are properly contracted with the polarization vector of the decaying particle to give a Lorentz scalar. This requires in some cases to include the tensor ε αβγδ p δ R via where ε αβγδ is the Levi-Civita symbol and j refers to the arguments of X defined in eqs. (4.22) and (4.23). Its antisymmetric nature ensures the correct parity transformation behavior of the amplitude. The spin factor for a whole decay chain, for example R → (R 1 → AB) (R 2 → CD), is obtained by combining the two-body terms and performing a sum over all unobservable, intermediary spin projections The main difference to the formalism used in ref. [13] is the inclusion of additional projection operators, i.e. P (S AB ) (R) and the one intrinsic to L (L AB ) (R), which ensure pure spin and angular momentum tensors. The spin factors for all decay topologies considered in this analysis are explicitly given in appendix B.

Measurement quantities
Here, we define all quantities derived from the amplitude model that are of physical importance. In order to provide implementation-independent measurements in addition to the complex coefficients a i , we define two quantities. Firstly, the fit fractions which are a measure of the relative strength between the different transitions. Secondly, the interference fractions are given by which measures the interference effects between amplitude pairs. Constructive interference leads to I ij > 0, while destructive interference leads to I ij < 0. Note that The global fractional CP -even content is defined as, is the decay amplitude for a D meson in a CP -even / CP -odd state. The parameter F + , can be determined from an amplitude model (eq. (4.28)) or by using model-independent methods [51]; the consistency of the two techniques provides a useful cross-check of the amplitude model. The fractional CP -even content also provides useful input to the determination of the CKM phase γ (φ 3 ) in B ± → DK ± and related decays. Additionally, knowledge of F + for all D decay final states can be used to determine the net CP -content of the D meson system, which is related to the charm-mixing parameter y D [52]. Finally, measurements of direct CP violation will also be reported. For this purpose, the amplitude coefficients are expressed in terms of a CP -conserving (c i ) and a CP -violating For ∆c i = 0 there is no CP violation between the corresponding D 0 and D 0 intermediate state amplitudes. Note that the CP -violating parameters are included only for distinct weak decay processes as the strong interaction is assumed to be CP -conserving such that e.g. the amplitudes for the processes D 0 → π − a 1 (1260) + → π + ρ(770) 0 and D 0 → π − [a 1 (1260) + → π + σ] share a common ∆c i , while having different CP -conserving parameters. As we do not measure the time distribution, we have no sensitivity to the overall phase difference between D 0 and D 0 and thus, the phase difference between A D 0 (x) and A D 0 (x) is fixed to null. From these separate amplitudes, the direct CP violation in each amplitude is simply calculated from the fit coefficients as In principle, the global direct CP asymmetry can be calculated from however to avoid an unnecessary systematic uncertainty arising from the amplitude model, this will instead be determined from an asymmetry in the integrated decay rates, composed of the number of signal candidates tagged as D 0 (D 0 ) mesons, N D 0 (N D 0 ). For the CLEO-c data, the signal tagging efficiency ratio, ε Tag ε Tag = 0.9899 ± 0.0015, (4.33) has been determined from an average over the D → Kπ, Kππ 0 and Kπππ efficiencies given in ref. [53]. No asymmetry in pion identification is found in the preceding CLEO data samples and thus the tagging efficiency ratio is set to unity with an uncertainty of 1.5% [54].

Likelihood fit
Due to flavor tagging, there are two independent data sets available; events which can be described by the amplitudes A D 0 (x) and A D 0 (x), respectively. In general, the signal PDF for events tagged as D 0 → h + h − π + π − is given by where Sig (x) is the phase-space efficiency and w is the wrong tag fraction as defined in section 3. In the case of no CP violation, the integrals over the D 0 and D 0 amplitudes will be equal. For the CP -tagged data sets used in the D 0 → K + K − π + π − analysis, the signal PDFs are given in ref. [13]. We do not account for effects of neutral charm meson oscillations, as we expect these to be negligible in these analyses. Note that the efficiency in the numerator appears as an additive constant in the log L that does not depend on any fit parameters such that it can be ignored. However, the efficiency function still enters via the normalization integrals. These normalization terms are determined numerically by a MC integration technique. For this purpose, we use simulated events generated according to a preliminary model, pass them through the full detector simulation and apply the same selection criteria as for data in order to perform the MC integrals. For example, the first integral in eq. (5.1) can be approximated as where A D 0 labels the preliminary amplitude model and x k is the k-th MC event. As a result, the efficiency can be included in the amplitude fit without explicitly modeling it. For D 0 → π + π − π + π − , we use a sample of N MC = 600000 MC events to ensure that the uncertainty on the integral is less than 0.5%. For D 0 → K + K − π + π − , we use samples of N MC ≈ 900000 events each, produced under each of the CLEO III and CLEO-c detector conditions. MC representing the CLEO II.V detector conditions is simulated from CLEO III MC via the reweighting process discussed in ref. [13]. The uncertainty on the integral for each D 0 → K + K − π + π − MC sample is less than 0.5%.
The background PDF, is determined in section 5.1 from sideband data. Note that because of the integration method, the background parameters only have meaning relative to the signal efficiency. The event likelihood is constructed from the signal PDF and the background PDF, where f Sig is the signal fraction as determined in section 3.1 and θ is the set of fit parameters.

Background model
Background events arise from randomly combined particles from various processes such as other D decays or continuum which, by chance, fulfill all required selection criteria. Some of them may even contain resonances that do not arise from the signal D 0 decay. The chosen background PDF for the D 0 → π + π − π + π − mode includes Breit-Wigner (BW) contributions from the resonances σ, ρ(770) 0 , f 0 (980) and two ad-hoc scalar resonances (S 0 1 , S − 2 ) with free masses and widths. They are added incoherently on top of two non-resonant components. In addition, several exponential and polynomial functions are included to allow for more flexibility. The background function is explicitly given by where, The real parameters b i , α i , c i , d i and e i are extracted from a fit to the sideband samples defined in section 3.1.

Signal model construction
The light meson spectrum comprises multiple resonances which are expected to contribute to D 0 → h + h − π + π − decays as intermediate states. Apart from clear contributions coming from resonances such as a 1 (1260) → ρ(770) 0 π, φ(1020) and K * (892) 0 , the remaining structure is impossible to infer due to the cornucopia of broad, overlapping and interfering resonances within the phase space boundary. The complete list of considered amplitudes can be found in appendix C.
To build the amplitude model, one could successively add amplitudes on top of one another until a reasonable agreement between data and fit was achieved. However, this stepwise approach is not particularly suitable for amplitude analyses as discussed in ref. [55]. Instead, we include the whole pool of amplitudes in the first instance and use the Least Absolute Shrinkage and Selection Operator [55,56] (LASSO) approach to limit the model complexity. In this method, the event likelihood is extended by a penalty term which shrinks the amplitude coefficients towards zero. The amount of shrinkage is controlled by the parameter λ, to be tuned on data. Higher values for λ encourage sparse models, i.e. models with only a few non-zero amplitude coefficients. The optimal value for λ is found by minimizing the Bayesian information criteria [57] (BIC), where N Sig is the number of signal events and r is the number of amplitudes with a decay fraction above a certain threshold. In this way, the optimal λ balances the fit quality (−2 log L) against the model complexity. The LASSO penalty term is only used to select the model. Afterwards, this term must be discarded in the final amplitude fit with the selected model, otherwise the parameter uncertainties would be biased. The implementation of the LASSO procedure differs between the D 0 → h + h − π + π − analyses. For D 0 → π + π − π + π − decays, the set of amplitudes is selected using the optimal value of λ = 28, and is henceforth called the LASSO model; figure 5(a) shows the distribution of BIC values obtained by scanning over λ where we choose the decay fraction threshold to be 0.5%. It is important to note that there are certain groups of amplitudes with the same angular distribution that are prone to produce artificially high interference effects. Amongst them are the di-scalar amplitudes: and D → f 0 (1370) f 0 (1370) as well as the di-vector amplitudes: D → (π π) P (π π) P , D → (π π) P ρ(1450) 0 and D → ρ(1450) 0 ρ(1450) 0 . In these cases, only one amplitude of the group is included at a time and the model selection is performed for each choice. It was further observed that the inclusion of the D → π[π(1300) → π(π π) P ] amplitude leads to a D → ρ(770) 0 ρ(770) 0 D-wave fraction much larger than the S-wave fraction with a large destructive interference. As we consider this as unphysical we do not include it in our default approach but in an alternative model presented in appendix D. In addition, we repeated the model selection procedure under multiple different conditions:  From that we obtained a set of alternative models shown in appendix D.
Due to the vast number of potential amplitude components and computational limits imposed by the consideration of multiple data samples in the D 0 → K + K − π + π − analysis, a staged LASSO method using only the flavor-tagged data, representing over 90% of the available statistics, is employed. The approach taken is based on the assumption that the signal decay proceeds primarily by doubly resonant decays, i.e. cascade and quasi-two-body decays, rather than decay amplitudes with non-resonant components. In Stage 1, only doubly resonant decays along with the simplest non-resonant component Figure 5(b) shows a plot of the complexity factor λ, against the resulting BIC values. We found that the fit cannot distinguish between amplitudes with K * (1680) + → K * (892) 0 π + and K * (1410) + → K * (892) 0 π + , which both peak outside the kinematic range of the D decay's phase space. We therefore only include K * (1680) + → K * (892) 0 π + in our nominal model. An alternative fit with the K * (1410) + , which has marginally worse fit quality is presented in table 20.
In Stage 2, the LASSO procedure is again performed with the components selected by Stage 1 and all single-resonant components. It should be noted in the case of cascade decays that if LASSO picked an amplitude component but not its conjugate decay in the -19 -

JHEP05(2017)143
first stage, the conjugate is also considered again in this stage. Once more, the interplay between D → SS amplitudes leads to very large interference terms, and thus f 0 (980) (π + π − ) S and f 0 (980) (K + K − ) S components are considered as a replacement for the non-resonant (K + K − ) S (π + π − ) S component in an alternative model. The final fit merges the components chosen in Stage 1 and Stage 2 and includes the CP -tagged data. Within this set of amplitudes, 6 are considered insignificant relative to their error and removed from the fit with no significant impact on fit quality.
6 D 0 → π + π − π + π − amplitude analysis results Table 3 lists the real and imaginary part of the complex amplitude coefficients a i , obtained by fitting the LASSO model to the data, along with the corresponding fit fractions. The letters in square brackets refer to the relative orbital angular momentum of the decay products. If no angular momentum is specified, the lowest angular momentum state consistent with angular momentum conservation and, where appropriate, parity conservation is used. The interference fractions are given in appendix E. Figure 6 shows the distributions of selected phase space observables, which demonstrate reasonable agreement between data and the fit model. We also project into the transversity basis to demonstrate good description of the overall angular structure in figure 7: the acoplanarity angle χ, is the angle between the two decay planes formed by the π + π − combination with minimum invariant mass, min[m(π + π − )], and the remaining π + π − combination in the D rest frame; boosting into the rest frames of the two-body systems defining these decay planes, the two helicity variables are defined as the cosine of the angle, θ, of each π + momentum with the D flight direction.

Amplitude model fit results
In order to quantify the quality of the fit in the five-dimensional phase space, a χ 2 value is determined by binning the data; where N b is the number of data events in a given bin, N exp b is the event count predicted by the fitted PDF and N bins is the number of bins. An adaptive binning used in ref. [13] is used to ensure sufficient statistics in each bin for a robust χ 2 calculation. At least 25 events per bin are required. The number of degrees of freedom ν, in an unbinned fit is bounded by N bins − 1 and (N bins − 1) − N par , where N par is the number of free fit parameters. We use the χ 2 value divided by ν = (N bins − 1) − N par as a conservative estimate. For the LASSO model, this amounts to χ 2 /ν = 1.40 with ν = 221 and N par = 34, indicating a decent fit quality.
In addition to the best five models as determined by the LASSO procedure, a further four alternative models are studied and presented in table 18. These comprise an "Extended" model whereby all conjugate partners of non-self-conjugate intermediate states chosen by the LASSO procedure are included. Two involving the removal of the π(1300) -20 -

JHEP05(2017)143
Decay channel 19.80 ± 3.54 ± 5.90     Figure 6. Invariant mass distributions of D 0 → π + π − π + π − signal candidates (points with error bars) and fit projections (black solid line). The signal component is shown in blue (dashed), the background component in red (dashed) and the wrongly tagged contribution in green (dashed). While the m 2 (π + π − ) includes all four possible π + π − combinations, the min[m 2 (π + π − )] (max[m 2 (π + π − )]) distribution includes the two π + π − combinations with the lowest (highest) invariant mass. The min[m 2 (π + π ± π − )] (max[m 2 (π + π ± π − )]) distribution includes the π + π ± π − combination with the lowest (highest) invariant mass. The effect of the K 0 S veto can clearly be seen in the top left projection. and a 1 (1640) resonances are described in the next section, while another based on the FOCUS model [12] is also considered. From this sample of alternative models, except the one based on the FOCUS model due to its poor fit quality, a model-dependent error on the fit fractions and the resonance parameters is derived from the variance. If one of the nominal amplitudes is not included in an alternative model, the corresponding fraction is set to zero. The dominant contribution is the a 1 (1260) resonance in the decay modes a 1 (1260) → ρ(770) 0 π and a 1 (1260) → σπ followed by the quasi-two-body decays D → σf 0 (1370) and D → ρ(770) 0 ρ(770) 0 . We find that the decay D 0 → a 1 (1260) + π − dominates over D 0 → a 1 (1260) − π + , which is similar to the pattern observed in the B sector, where B 0 → a 1 (1260) + π − is preferred over B 0 → a 1 (1260) − π + [59, 60].
It is important to note that even though the a 1 (1640) and the π(1300) resonances are selected by the model building, satisfactory fit results can also be obtained without them. The LASSO models obtained when explicitly excluding the a 1 (1640) and the π(1300) resonance from the pool of amplitudes are given in appendix D. These models are used to generate many pseudo-data sets according to the "no-a 1 (1640)" or "no-π(1300)" hypotheses denoted as H 0 . The pseudo-data is then fitted with H 0 and the alternative hypotheses, e.g. a 1 (1640) hypothesis H 1 , in order to predict the distributions of the log-likelihood differences ∆(−2 log L) = 2 log(L(H 1 )/L(H 0 )) under the H 0 hypotheses. We use a Gaussian function to parameterize the ∆(−2 log L) distributions. By integrating the tails of the Gaussians above the ∆(−2 log L) value observed on the real data, the H 0 hypotheses can be excluded in favor of the a 1 (1640) and π(1300) alternate hypotheses at the 2.4σ and 6.1σ levels, respectively.
Since the a 1 (1640) + resonance is not yet well established, we verify its resonant phase motion in a quasi-model-independent way as pioneered in ref. [68]. For this purpose, the Breit-Wigner lineshape is replaced by a complex-valued cubic spline. The interpolated cubic spline has to pass through independent complex knots spaced in the m 2 (π + π + π − ) region around the nominal mass. The position of the knots is chosen ad-hoc. We verified on simulated experiments that with this choice a Breit-Wigner lineshape can be properly reproduced, given there is a real resonance. The fitted magnitudes and phases of the knots are shown in figure 8, where the expectations from a Breit-Wigner shape with the mass and width from the nominal fit are superimposed taking only the statistical uncertainties on the mass and width into account. The interpolated spline generally reproduces the features of the Breit-Wigner parametrization. In particular, the resulting   Argand diagram shows a circular, counter-clockwise trajectory which is the expected behavior of a resonance. Note that the high-mass tail of the a 1 (1640) is outside of the phase space boundary such that it is not possible to investigate the full phase motion. Similar quasi-model-independent studies are performed for the a 1 (1260) and π(1300) resonances as shown in figures 9 and 10, respectively. Since the investigated resonances are all very broad, the quasi-model-independent lineshapes can absorb statistical fluctuations in the data, especially near the phase space boundaries. Therefore, the agreement with the Breit-Wigner expectation in all cases indicates that it is qualitatively reasonable that these resonances are indeed real features of the data.

Global CP content measurement
The fractional CP -even content, F 4π + , is determined from the integral in eq. (4.28), using the nominal model for A 4π D 0 and A 4π D 0 assuming no direct CP violation in the D meson decay. The uncertainty on F 4π + is calculated from pseudo-experiments by randomly varying the free parameters of the amplitude fit within their measured statistical and systematic uncertainties. For each variation, F 4π + is redetermined, and the square root of the sample variance of these values is taken as the uncertainty. An additional systematic uncertainty is assigned by computing F 4π + for each of the alternative amplitude models. The standard -26 -
deviation of these values is taken as the additional model uncertainty. The obtained result,

Search for direct CP violation
A search for CP violation is performed by fitting the LASSO model to the flavor-tagged D 0 and D 0 samples. In contrast to our default fit described in section 5, we now allow the amplitude coefficients for D 0 → π + π − π + π − and D 0 → π − π + π − π + decays to differ, as described in section 4.3.
The masses and widths of the resonances are fixed to the values obtained in the nominal fit. Possible additional biases due to this assumption are included in the systematic uncertainties which are otherwise determined as described in section 8. Table 6 compares the resulting fit fractions for the D 0 and D 0 decays. The sensitivity to A i CP is at the level of 4% to 22% depending on the decay mode. No significant CP violation is observed for any of the amplitudes. Also, the integrated CP asymmetry over phase space is found to be which is consistent with CP conservation. Due to the cancellation of systematic uncertainties in asymmetry-like quantities, the only remaining source considered for the global CP asymmetry is the tagging efficiency ratio, which is set to unity for this purpose. This nominal value of A 4π CP is consistent with that which can be found from the amplitude model via eq. 7 D 0 → K + K − π + π − amplitude analysis results Table 7 lists the real and imaginary part of the complex amplitude coefficients a i , along with the corresponding fit fractions. The interference fractions are given in appendix E. Figures 11 and 12 show the distributions of selected phase space observables, which demonstrate reasonable agreement between data and the fit model. For the flavor-tagged data only, we also project into the transversity basis to demonstrate good description of the overall angular structure in figure 13: the acoplanarity angle χ, is the angle between the two decay planes formed by the K + K − combination and the π + π − combination in the D rest frame; boosting into the rest frames of the two-body systems defining these decay planes, the two helicity variables are defined as the cosine of the angle, θ K + , of the K + momentum with the D flight direction, and the cosine of the angle, θ π + , of the π + momentum with the D flight direction. In contrast to the treatment of the a 1 (1260) and π(1300) substructure in the D 0 → π + π − π + π − analysis, we do not enforce the same amplitude substructure for the K(1270) + , K(1400) + , K(1680) + decays as for K(1270) − , K(1400) − , K(1680) − ; this choice has historical reasons. It is re-assuring to see that the results we obtain without these constraints are consistent with what one would expect if such constraints had been applied (cf. model A in table 20). For the LASSO model, the χ 2 /ν is 1.5 with ν = 116, where the effective number of degrees of freedom is determined with a pseudo-experiment technique. Its value is chosen to be the one that best converts the distribution of χ 2 values for each experiment into the standard uniform distribution. This method differs from that used in D 0 → π + π − π + π − as the relatively small size of the data sample here would otherwise result in negative degrees of freedom.

Amplitude model fit results
Four alternate models are presented in appendix D: (A) a model that requires the use of conjugate amplitudes for all present non-selfconjugate decays (B) replacing K * (1680) + → K * (892) 0 π + with the K * (1410) + → K * (892) 0 π + amplitude (C) replacing the flat non-resonant term with the f 0 (980) (π + π − ) S and f 0 (980) (K + K − ) S amplitudes (D) the model previously reported in ref. [13] The results between models are broadly consistent where the largest individual fit fraction corresponds to the D 0 → φ(1020) ρ(770) 0 amplitude. We found that we cannot distinguish between the K * (1680) meson in our default model and the K * (1410) meson trialled in alternative model B. Both of these components peak outside the kinematically allowed range.
Relative to the previous analysis of the same data set [13], the most notable apparent difference in our default model is the fit fraction of the φ(1020) ρ(770) 0 S-wave, which was 38.3% in ref. [13], but only 28.1% in our current analysis. This is because of our modified description of the V V D-wave. In ref. [13], the component labeled as D-wave is -28 -

JHEP05(2017)143
Decay channel 6.36 ± 1.24 ± 3.41 -6.86 ± 1.37 ± 3.38 34.93 ± 3.74 ± 8.39 5.66 ± 4.50 ± 6.92 12.4 ± 2.6 ± 3.9 ± 5.0 43.52 ± 4.59 ± 9.73 -1.20 ± 5.20 ± 5.27 29.62 ± 7.39 ± 12.19 -54.93 ± 5.65 ± 11.92 0.42 ± 0.05 ± 0.10 0.49 ± 0.05 ± 0.09 11.1 ± 1.2 ± 2.1 ± 0.7 Sum 106.9 ± 4.5 ± 6.9 ± 6.1 Table 7. Real and imaginary part of the complex amplitude coefficients and fit fraction of each component of the D 0 model. For the fit coefficients, the first quoted uncertainty is statistical, while the second arises from systematic sources. The third uncertainty in the fit fraction arises from the alternative models considered.  Figure 11. Invariant 2-body mass distributions of D 0 → K + K − π + π − signal candidates shown as points with error bars. The overall fit projection is shown in black, the signal in blue and the background in red. The effect of the K 0 S veto can clearly be seen in the bottom right projection. a superposition of D and S waves, a choice which was motivated by the convention used in four-body amplitude analyses at the time. This led to a large interference between the components labeled as S wave and D wave of -15.7%. In this analysis, as we parametrize a pure D-wave, we find an interference fraction between the φ(1020) ρ(770) 0 S-and D-waves of -3.7%. Taking these interference fractions into account, the combined φ(1020) ρ(770) 0 Sand D-wave fraction of 26% is therefore consistent between both analyses. In contrast to ref. [13], we also find a small, but significant φ(1020) ρ(770) 0 P -wave component. Another difference in the two-resonance topology is in the K * (892) 0K * (892) 0 mode, where our results indicate a significant P -and D-wave contribution, while in ref. [13], only an S- wave contribution was observed. Note though, that model 6 in ref. [13] has a P -wave in the K * (892) (Kπ) P decay of a similar size as our K * (892) 0K * (892) 0 P -wave. The largest differences in our results are, as might be expected, in the cascade topology, because of the significant changes we implemented to improve the description of the lineshapes of resonance decays to three-body final states. We find that the process D 0 → K * * + K − , where K * * represents any excited kaon, dominates over D 0 → K * * − K + , analogous to the dominance of D 0 → a 1 (1260) + π − over D 0 → a 1 (1260) − π + decays. In ref. [13], this was only the case for the K(1270) → K * (890) π amplitude. We also observe a significant K(1270) → K * (1430) π component in agreement with ref.
[14] but not with ref. [13]. The description of this type of decay chain, with a daughter whose mean mass is outside the kinematically allowed region, benefits particularly from our improved lineshapes. As in ref.

Global CP content measurement
Following the same approach as for D 0 → π + π − π + π − decays, for the fractional CP -even content we obtain

Search for direct CP violation
Following the same approach as for D 0 → π + π − π + π − decays, we measure the direct CP violating parameters given in All measurements are consistent with CP conservation.

Systematic uncertainties
There are three main sources of systematic uncertainties on the fit parameters to be considered; an intrinsic fit bias, as well as experimental and model-dependent uncertainties. For each four-body decay, the fit bias itself is determined from a large ensemble of MC pseudoexperiments generated from the nominal LASSO model. The mean difference between the generated and fitted parameters are taken as a systematic uncertainty. The experimental systematic uncertainties occur due to imperfect knowledge of the yield of background events and their distribution in phase space, the wrong tag probability, and various effects on the efficiency variation over phase space. To estimate the systematic uncertainty related to the background shape that was fixed from sideband, the -32 -
amplitude fit is repeated where the background parameters are allowed to vary within their statistical uncertainties. In addition, several alternative background PDFs are tested whereby each background contribution is replaced, one at a time, by a flat, non-resonant model. The largest deviations from the nominal values are assigned as systematic uncertainties. The uncertainty due to both the signal fraction and the wrong tag probabilities in the flavor-tagged samples are estimated by repeating the fit and allowing them to vary under Gaussian constraints. The signal fraction uncertainty for the CP -tagged sample is determined by fixing the fraction to unity and repeating the fit. Various assumptions made on the acceptance in the fit model are also considered. As the acceptance comes from MC, we account for differences between data and MC arising from tracking and particle identification as a function of momentum of the daughter particles. Using correction factors obtained from independent internal CLEO studies, the MC is reweighted separately for each effect and the fit to data repeated. While detector resolution can be safely ignored in D 0 → π + π − π + π − decays, neglecting the effect of finite momentum resolution on the φ(1020) resonance in D 0 → K + K − π + π − decays may lead to a bias. To counter this, a large number of pseudo-experiments were generated by distributing MC events that have passed full selection and weighted by the LASSO model found from data. Each experiment is then fit with the signal model where the mean difference between the generated and fitted parameters are assigned as systematic uncertainties. Finally, the integration error due to the limited size of the MC sample is of the order of 0.5%, so it is neglected as a source of systematic uncertainty.
Model-dependent uncertainties arise from fixed lineshape parameters and the effects of interference from Cabbibo-suppressed decays on the tag-side in the CLEO-c flavortagged data samples. The uncertainties due to fixed masses and widths of resonances are evaluated by varying them one-by-one within their quoted errors. In our nominal fit, the Blatt-Weisskopf radial parameter is set to r BW = 1.5 (GeV/c) −1 . As a systematic check, we set the radial parameter to zero. For the calculation of the energy-dependent widths, -33 -

JHEP05(2017)143
the partial widths into the πππ channel are obtained using an iterative procedure described in section 4.1. The systematic error of this approach is estimated by repeating the fit using the iteration previous to the final. In some cases, the energy-dependent width relies on external measurements of intermediate branching fractions. In D 0 → π + π − π + π − , their impact is studied by recalculating the width considering only decays into the πππ (ππ) final state for three-body (two-body) resonances. For D 0 → K + K − π + π − , the energydependent widths of the three-body resonances are recalculated assuming a flat phase space distribution. Similarly, the energy-dependent mass of the a 1 (1260) resonance is approximated by a constant and the resulting shifts of the fit parameters are assigned as systematic errors.
The systematic uncertainty related to interference from the tag-side arising between the CKM-favored c → s and CKM-suppressed c → d amplitudes in the final states used for flavor-tagging is accounted for by using an alternative signal PDF at the cost of two additional fit parameters as described in ref. [13].

Conclusion
The first amplitude analysis of flavor-tagged D 0 → π + π − π + π − decays has been presented based on CLEO-c data. Due to the large amount of possible intermediate resonance components, a model-building procedure has been applied which balances the fit quality against the number of free fit parameters. The prominent contribution is found to be the a 1 (1260) resonance in the decay modes a 1 (1260) → ρ(770) 0 π and a 1 (1260) → σ π. Along with the a 1 (1260), further cascade decays involving the resonances π(1300) and a 1 (1640) are also seen. The masses and widths of these resonances are determined using an advanced lineshape parametrization taking into account the resonant three-pion substructure. The resonant phase motion of these states has been verified by means of a quasi-model-independent study. In addition to these cascade topologies, there is a significant contribution from the quasi-two-body decays D 0 → ρ(770) 0 ρ(770) 0 and D 0 → σ f 0 (1370). The CP -even fraction of the decay D 0 → π + π − π + π − as predicted by the amplitude model is consistent with a previous model-independent study. The amplitude model has also been used to search for CP violation in D 0 → π + π − π + π − decays, where no CP violation among the amplitudes is observed within the given precision of a few percent.
Moreover, the amplitude analysis of D → K + K − π + π − decays performed by CLEO [13] has been revisited by applying the significantly improved formalism presented in this paper, using decays obtained from CLEO II.V, CLEO III, and CLEO-c data. The largest components are the processes D 0 → φ(1020) ρ(770) 0 , D 0 → K 1 (1270) + K − and D 0 → K(1400) + K − , which together account for over half of the D 0 → K + K − π + π − decay rate. The fractional CP -even content is measured for the first time and a search for CP asymmetries in the amplitude components yields no evidence for CP violation.
In addition to shedding light on the dynamics of D 0 → h + h − π + π − decays, these results are expected to provide important input for a determination of the CP -violating phase γ (φ 3 ) in B − → DK − decays.  Table 9. Systematic uncertainties on the fit parameters of our nominal D 0  Table 11. Systematic uncertainties on the fit parameters of our nominal D 0  A Energy-dependent widths

B Spin amplitudes
The spin factors used for D → P 1 P 2 P 3 P 4 decays are given in table 13. To fix our phase convention, we give the exact matching of the particles P 1 , P 2 , P 3 , and P 4 in the spin factor definition to the final state particles in specific decay chains in tables 14 and 15.

Number
Decay chain Spin amplitude Table 13. Spin factors for all topologies considered in this analysis. In the decay chains, S, P , V , A, T and P T stand for scalar, pseudoscalar, vector, axial vector, tensor and pseudotensor, respectively. If no angular momentum is specified, the lowest angular momentum state compatible with angular momentum conservation and, where appropriate, parity conservation, is used.
Decay channel Spin factor number P 1 Table 15. Spin factors used for the decay chains included in the D 0 → K + K − π + π − LASSO model, including the particle numbering scheme. The second column refers to the spin factors as numbered in table 13, and the particles P 1 , P 2 , P 3 , and P 4 refer to those defined in table 13.

JHEP05(2017)143 C Considered decay chains
The various decay channels considered in the model building are listed in tables 16 and 17.
Supplemental material. We provide a collection of C macros to reproduce all energydependent masses and widths described in section 4.1. These are intended to be parsed by the ROOT software and have names indicating which energy-dependent quantity and resonance they correspond to.
Two additional text files containing the statistical correlation matrices of the nominal results for D 0 → π + π − π + π − and D 0 → K + K − π + π − are provided. Their filenames are Correlations4pi.txt and CorrelationsKKpipi.txt, respectively. The format of each file is as follows. Firstly, each free parameter is assigned a numerical identifier. Following this, the lower diagonal correlation matrix is given for these indices.