Search for $CP$ violation through an amplitude analysis of $D^0 \rightarrow K^+ K^- \pi^+ \pi^-$ decays

A search for $CP$ violation in the Cabibbo-suppressed $D^0 \rightarrow K^+ K^- \pi^+ \pi^-$ decay mode is performed using an amplitude analysis. The measurement uses a sample of $pp$ collisions recorded by the LHCb experiment during 2011 and 2012, corresponding to an integrated luminosity of 3.0 fb$^{-1}$. The $D^0$ mesons are reconstructed from semileptonic $b$-hadron decays into $D^0\mu^- X$ final states. The selected sample contains more than 160000 signal decays, allowing the most precise amplitude modelling of this $D^0$ decay to date. The obtained amplitude model is used to perform the search for $CP$ violation. The result is compatible with $CP$ symmetry, with a sensitivity ranging from 1% to 15% depending on the amplitude considered.


Introduction
The asymmetry between matter and antimatter in the universe is one of the most compelling questions that awaits explanation in particle physics.One of the three conditions required to explain this asymmetry is charge-parity (CP ) violation [1].The Standard Model (SM) of particle physics allows this violation to arise.This is, however, not sufficient to describe the cosmological observations when taking into account the estimated lifetime of the universe [2,3].
A promising area to search for CP violation beyond the SM is the decay of charm hadrons.In Cabibbo-suppressed charm-hadron decays, the interference between the c → q (q = d, s) tree diagrams with the c → u loop diagrams is the only potential source of CP asymmetry in the SM, which is predicted to be smaller than 0.1% [4,5].These decays are sensitive to potential new contributions in strong penguin and chromomagnetic dipole operators that could produce effects of up to 1% [6].Multibody decays of charm hadrons have peculiarities that make them particularly interesting for these studies.They have a rich resonant structure, and the variation of the strong phases over the decay phase space may provide regions with enhanced sensitivity to CP violation.
This paper presents a search for CP violation in the individual amplitudes of the decay1 D 0 → K + K − π + π − , produced in semileptonic decays of b hadrons into D 0 µ − X, where X represents any combination of undetected particles.The charge of the µ − (µ + ) lepton is used to identify D 0 (D 0 ) mesons.The D 0 µ − system is reconstructed in a sample of pp collisions collected by the LHCb experiment in 2011 and 2012 at centre-of-mass energies of 7 and 8 TeV, corresponding to integrated luminosities of 1.0 fb −1 and 2.0 fb −1 , respectively.A CP -averaged decay model is developed through a full amplitude analysis in the five-dimensional phase space of the four-body D 0 decay, where the D 0 and D 0 samples are merged.This model is then used to search for CP violation in a simultaneous fit to the distinct D 0 and D 0 samples.
Such an amplitude analysis has already been performed by the CLEO experiment [7] and has been recently updated with the CLEO legacy data [8].The analysis was based on 3000 signal decays and no CP violation was observed.This paper presents an analysis using a data sample more than 50 times larger.Searches for CP violation in this decay have also been performed by the BaBar [9], LHCb [10] and Belle [11] experiments, which are all compatible with CP symmetry.
In Sec. 2, a description of the LHCb detector and the simulation is presented.The event selection procedure is described in Sec. 3. The formalism used for the amplitude analysis is presented in Sec. 4. The systematic uncertainties are discussed in Sec. 5.The resulting CP -averaged model is presented in Sec. 6.The CP -violation results are summarised in Sec.7 and conclusions are discussed in Sec. 8.

Detector and simulation
The LHCb detector [12,13] 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 siliconstrip 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 of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c.The polarity of the dipole magnet is reversed periodically throughout data-taking.The configuration with the magnetic field vertically upwards (downwards) bends positively (negatively) charged particles in the horizontal plane towards the centre of the LHC.The minimum distance of a track to a primary vertex (PV), the impact parameter (IP), 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 detectors.Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic and a hadronic calorimeter.Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.The online event selection is performed by a trigger consisting of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction.
In the simulation, pp collisions are generated using Pythia [14] with a specific LHCb configuration [15].Decays of hadronic particles are described by EvtGen [16], in which final-state radiation is generated using Photos [17].The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [18] as described in Ref. [19].

Event selection
At the hardware trigger stage, events are required to have a muon with high p T or a hadron with high transverse energy in the calorimeters.For the muon, the p T threshold is 1.48 GeV/c in 2011 and 1.76 GeV/c in 2012.For hadrons, the transverse energy threshold is 3.5 GeV and 3.62 GeV, respectively.The software trigger requires a two-, three-or four-track secondary vertex with a significant displacement from any pp interaction vertex.At least one charged particle must have a transverse momentum p T > 1.5 GeV/c and be inconsistent with originating from a PV.A multivariate algorithm [20] is used for the identification of secondary vertices consistent with the decay of a b hadron.
In the offline selection, two kaons and two pions are combined to form a D 0 candidate, which is itself combined with a muon to form a b-hadron candidate.These candidates are required to pass preselection criteria, listed below.The kaons and pions are required to have a momentum larger than 2 GeV/c and a transverse momentum larger than 0.3 GeV/c.The muon is required to have a momentum larger than 3 GeV/c and a transverse momentum larger than 1.2 GeV/c.Requirements on the track quality as well as correct particle identification are also made.The D 0 candidate is required to decay downstream of the b-hadron decay vertex and its mass is required to be in the range [1805,1925] MeV/c 2 .The mass of the b-hadron candidate is required to be in the range [2500,6000] MeV/c 2 .Furthermore, events are kept only if the positive trigger decision is due to the signal candidate.
A boosted decision tree (BDT) [21,22] is used to separate signal from background.The description of the signal is taken from simulation while that of the background is taken from the D 0 -mass sidebands in data.The two datasets are split into two independent subsamples for the training and the testing of the classifier.Many kinematic and isolation variables are tested, the set of variables being gradually reduced until the performance of the BDT starts to decrease significantly.The final set of variables used in this analysis contains the mass of the b hadron corrected for the missing particles [23], the flight distance of the D 0 candidate, the vertex-fit χ 2 of the decay vertex of the D 0 candidate, particle identification information of the four daughters of the D 0 candidate, the probability of the tracks being reconstructed from random hits [13], and isolation variables based on the underlying event.One isolation variable is obtained by considering all tracks inside a cone around the D 0 direction and computing the p T asymmetry with respect to the D 0 meson.Two additional isolation variables are obtained by adding tracks to the D 0 vertex; the first represents the invariant mass of the new combination when adding the track with the minimum vertex χ 2 difference and the second represents the minimum χ 2 difference when adding two tracks.The selection on the BDT outcome is chosen to maximise the figure of merit N s / √ N s + N b , where N s and N b are the number of signal and background events in the signal region (SR), where the SR is defined as ±2 standard deviations around the central value of the known D 0 mass [24].
A background component arises when the D 0 meson originates from a D * + resonance, which is itself coming from a b-hadron decay, and the slow pion from the D * + resonance is used as one of the pions in the D 0 decay.This misreconstruction generates a structure in the upper end of the K + K − π − invariant mass spectrum, which is removed with the requirement m(K is present in the SR.These decays have a different topology than D 0 → K + K − π + π − decays and are therefore vetoed by removing all candidates that have a π + π − invariant mass in the region [480.2,507.2]MeV/c 2 .Background contributions in which both a kaon is misidentified as a pion and a pion as a kaon are found to be negligible.About 1.7% of the events contain more than one D 0 candidate, in which case one candidate per event is randomly selected.The resulting D 0 mass distribution, shown in Fig. 1, is fitted with a double Gaussian signal function and an exponential background function.Only the N data = 196 648 events that are in the SR are kept.According to the fit, the selected sample contains 162 909 ± 516 signal decays in the SR with a purity of f s = 82.8± 0.3%.This selection provides a significant improvement with respect to the cut-based selection of a previous LHCb analysis of the same decay mode with the same data sample [10]: 5% more signal and 30% less background are retained.
In order to improve the resolution on the measured momenta of the D 0 decay products, the tracks and vertices are refitted under the constraint that the invariant mass of the four D 0 daughter particles be equal to the known D 0 mass [24].
The data sample contains both decays, tagged by the sign of the accompanying muon.For the rest of the analysis, the CP transformation is applied on the D 0 candidates, i.e. on the momentum vectors of the four daughter particles in the D 0 rest frame.This allows, in absence of CP violation, the D 0 and CP -transformed D 0 candidates to have identical distributions.

Amplitude analysis
The amplitude analysis consists of describing the D 0 decay chain as a coherent sum of amplitudes, each corresponding to a specific decay path (called "component" from here on) from the mother particle to the K + K − π + π − final state.These complex amplitudes may interfere.The main goal of the amplitude analysis is to identify the components that contribute to the decay.Each amplitude A k is multiplied by a complex coefficient c k , whose modulus |c k | and phase arg(c k ) are determined by means of an unbinned maximum likelihood fit (except for one of the complex coefficients, which is fixed to 1).

Likelihood fit
The likelihood function is defined as where a(x; c) and b(x) are the probability density functions (PDF) of the signal and background components normalized as a(x; c)d 5 x = b(x)d 5 x = 1.Here, c represents the fit parameters and x represents the five dimensions of the D 0 → K + K − π + π − four-body phase space.These five dimensions can be visualised with the five Cabibbo-Maksymowicz (CM) variables [25] (see Fig. 2): the invariant mass of the two-kaon system m(K + K − ); the invariant mass of the two-pion system m(π + π − ); the cosine of the helicity angle for the two-kaon system cos(θ K ), defined as the angle between the direction of the D 0 momentum and that of one of the kaons in the rest frame of the two kaons; the cosine of the helicity angle for the two-pion system cos(θ π ), defined similarly; and the angle φ between the plane defined by the directions of the two kaons and the plane defined by the directions of the two pions, in the D 0 rest frame.The signal PDF is written as where s (x) is the signal efficiency, S(x; c) is the signal model described in Sec.4.2 and R 4 (x) is the function representing the four-body phase space.Instead of explicit parametrisations of the signal efficiency and the four-body phase-space function, the fit uses a simulated signal sample that encodes the functions s (x) and R 4 (x).After the reconstruction and selection are applied, this sample is distributed according to the PDF a gen (x) = s (x)S gen (x)R 4 (x)/I gen , where S gen (x) is the signal model used in the generation of the sample and I gen = s (x)S gen (x)R 4 (x)d 5 x.This simulated sample contains approximately 10 million signal events, half of which are generated with a pure phase-space decay model (corresponding to a constant function S gen (x)) in order to cover the whole phase space and the other half according to the decay model published by the CLEO collaboration [7] in order to sufficiently populate the phase-space regions corresponding to the resonances.The normalisation integral of Eq. 2 is then estimated as where N sim is the number of events in the simulated sample.The background component is described using the same simulated sample.For each simulated event i at position x i in phase space, a weight w(x i ) is assigned so that the weighted simulated distribution matches the distribution of the background in five dimensions, as explained in Sec.4.3.These weights are obtained by evaluating w(x) = b(x)/a gen (x) and satisfy 1 N sim N sim i=1 w(x i ) = 1.In order to factorise out the signal efficiency and the phase-space function, the background PDF is rewritten as b(x) = s (x)B(x)R 4 (x)/I gen , where B(x) = w(x)S gen (x).
Using the above definitions, the likelihood is rewritten as where C = N data j=1 ( s (x j )R 4 (x j )/I gen ) is a constant independent of the fit parameters c.Therefore, this constant does not need to be computed in order to maximize L(c).The value B(x j ) of data event j is obtained as the average of the quantities B(x i ) over all simulated events with x i falling in the same phase-space bin as x j .For this purpose an adaptive binning is used, which splits five-dimensional rectangular bins along the CM variables until each bin contains between 5 and 9 simulated events, resulting in smaller bin volumes in higher density regions.

Signal description
The formalism2 chosen for this amplitude analysis is the so-called isobar model [26,27], which assumes that each amplitude can be built as a series of two-body decays.The two allowed patterns for D 0 → abcd decays, both involving two intermediate resonances r 1 and r 2 , are D 0 → r 1 r 2 followed by r 1 → ab and r 2 → cd, and D 0 → r 1 a followed by r 1 → r 2 b and r 2 → cd.In both cases the k th complex amplitude A k (x) is computed as the product of the lineshapes for the resonances r 1 and r 2 , the normalised Blatt-Weisskopf barrier factors [28], and a spin factor defined using the covariant formalism [29].
The resonances considered in this analysis are listed in Table 1.The default lineshape, used for most resonances, is the relativistic Breit-Wigner (RBW) function [30].For the a 1 (1260) + and the K 1 (1270) + resonances, a correction is applied to their mass-dependent width according to the formalism described in Ref. [31], in order to take into account the effect of intermediate resonances.In addition, an exponential form factor derived from Refs.[31,32] is used for these two three-body resonances instead of the Blatt-Weisskopf factors.Exceptions are the Flatté parametrisation [33] used for the a 0 (980) 0 resonance near the KK threshold and the Gounaris-Sakurai parametrisation [34] used for the ρ(770) 0 resonance.The aforementioned lineshapes describe accurately well-separated narrow resonances.In the case of broad overlapping resonances, the above descriptions may fail to account properly for interferences and the K-matrix formalism [35] is used instead.The ππ and the KK S-wave contributions (referred to as [π + π − ] L=0 and [K + K − ] L=0 in the following) are both described with the K-matrix formalism.They couple to five different channels (ππ, KK, ππππ, ηη and ηη ) with five different poles (f 0 (980), f 0 (1300), f 0 (1500), f 0 (1750) and f 0 (1200 − 1600)) and a non-resonant contribution, according to the parametrisation of Ref. [36], with parameters taken from Ref. [37].The Kπ S-wave contribution (referred to as [K + π − ] L=0 in the following) is also described with the K-matrix formalism.It couples to the Kπ and Kη channels in the isospin-state I = 3  2 , or only to the Kπ channel if I = 1  2 , with a single pole (K * 0 (1430)) and a non-resonant contribution.
The mass and the width of the K 1 (1270) + resonance are left floating in the final fit.The parametrisation of the a 1 (1260) + resonance is taken from Ref. [39] with mass and width set to 1195 MeV/c 2 and 422 MeV/c 2 .For all the other resonances, their masses and widths are fixed to the values of Ref. [24].The radii of the normalised Blatt-Weisskopf barrier factors used in the parametrisation of the D 0 , K * (892) 0 and K * (1680) 0 mesons are set to 1.21 GeV −1 , 1.13 GeV −1 and 1.93 GeV −1 , respectively.Each of these three values is obtained by maximising the likelihood of amplitude fits while fixing the values of the other two.The radius of the a 1 (1260) + resonance is fixed to 1.7 GeV −1 according to Ref. [39] and all the other radii are fixed to 1.5 GeV −1 .
The total signal function is then described by a coherent sum of N amplitudes, where c 1 = 1 and the other complex coefficients c k are defined relative to c 1 .The moduli and phases of these other complex coefficients are left floating in the fit.After the fit, in order to quantify the relative importance of each component, the fit fraction is computed.Interferences may lead to N k=1 F k = 1.

Background description
The sidebands of the D 0 -mass distribution are used to describe the background, assuming that the sum of the lower (m(K .92] GeV/c 2 ) sidebands gives a good description of the background in the SR.Although the background has a nontrivial distribution in the fivedimensional phase space, as it also contains resonances, this assumption is justified by the observation that the five-dimensional distributions of the candidates in the two sidebands are similar and that the contributions of the various resonances vary smoothly as a function of the Since the point in phase space of each candidate is recomputed under the D 0 -mass constraint, the phase-space boundaries are the same for candidates in the SR and in each of the sideband regions.However, as a side effect, the peaks of the resonances are slightly shifted to larger (smaller) values in the lower (upper) sideband.This effect is corrected for by weighting the sideband events to shift the resonance peaks back to their correct position while keeping their widths unchanged.The correction is applied for the φ(1020), K * (892) 0 and K * (892) 0 peaks visible in the K + K − , K − π + and K + π − invariant mass distributions, respectively.Since the K * (892) 0 peak is much less pronounced than the K * (892) 0 peak, the positions and shapes of these two peaks are constrained to be the same.No correction is applied for the small ρ(770) 0 contribution in the π + π − invariant mass of the sideband sample.
After these corrections, the sample of fully simulated signal events is reweighted to match the distribution of the sideband data.To obtain the weights w(x i ) needed for the amplitude fit (see Sec. 4.1), a multidimensional reweighting is performed simultaneously on 31 different projections of the five-dimensional phase space (10 invariant masses, 18 helicity angles, and 3 decay plane angles) using the package hep ml [40].

Model building
In order to build a CP -averaged decay model, a single model is fitted to the full data sample, containing both D 0 and CP -transformed D 0 candidates.This allows the construction of the signal model in a way that is blind to possible CP -violating effects.
The signal model is built iteratively.As a starting point, the model is made of the sum of the D 0 → φ(1020)(ρ − ω) 0 and D 0 → K * (892) 0 K * (892) 0 amplitudes, where in both cases three different values of the orbital angular momentum (L = 0, 1, 2) between the vector resonances from the D 0 decay are allowed.The ρ(770) 0 and ω(782) mesons, being close in mass, interfere heavily.They are therefore described as a quantum superposition of the two individual states, defined as the "(ρ − ω) 0 " state, with a free relative complex coefficient c.A long list of other possible amplitudes is defined from all possible combinations of the resonances listed in Table 1.At each iteration, the amplitude of this list that produces the largest decrease in the minimised value of −2 ln(L(c)) is added to the model.
As no CP violation is expected to arise in strong decays, the two charge states of the same three-body resonance are constrained to have the same decay substructure, implying that the two charged-conjugate states are always added together in the model building procedure.The K * (1410) and K * (1680) resonances give similar contributions and the fit cannot distinguish them.One of the two needs therefore to be removed.The K * (1680) components are chosen as the fit shows a slightly better χ 2 , but the K * (1410) components are considered in one of the alternative models used to estimate systematic uncertainties.
As iterations proceed, the sum of the fit fractions remains quite stable but diverges at some point.The procedure is then stopped and the model from the previous iteration is taken as the nominal one.

CP -violating observables
For measuring CP violation, the data is split in two samples according to the charge of the muon, one with D 0 decays, and another one with CP -transformed D 0 decays.This allows the use of the same model in a simultaneous fit to the two samples, where different parameter values are allowed for the two samples.Any significant difference between the fitted parameters on the two samples will signal CP violation.For each amplitude in the model, the fit is parametrised with the average modulus |c k |, modulus asymmetry A |c k | , average phase arg(c k ) and phase difference ∆ arg(c k ), defined as where |c k | and arg(c k ) are the polar coordinates (modulus and phase) of the complex fit coefficient multiplying the k th amplitude.The simultaneous fit minimises the sum of the two negative log-likelihoods for the D 0 and D 0 samples.Since no CP violation is expected in the strong decays of the three-body resonances, their modulus and phases are therefore simultaneously fitted to common values for the two samples.An additional information on CP violation in each amplitude can be obtained from the fit fractions.The asymmetry is considered, where are the fit fractions of the k th amplitude for the D 0 and the D 0 samples respectively.

Systematic uncertainties
Various systematic uncertainties are considered, most of them related to the fitting procedure and the model determination, and others to the interaction of particles with the detector.Among the effects that could directly influence the amplitude model are the determination of the selection efficiency, the background description, the signal fraction, the description of the resonance shapes, alternative components in the model, and any bias intrinsic to the fit procedure.Moreover, the flavour misidentification and the track reconstruction asymmetry may have an effect on the CP -violation measurements.
The size of these uncertainties is determined either by running pseudoexperiments (fit bias, background description, flavour misidentification, detection asymmetry, alternative models), by performing alternative fits after applying some modifications (background description, selection efficiency, resonance shapes), or by fitting multiple times the data to simulate statistical fluctuations using resampling (background description) or to take into account the uncertainties on the fixed parameters (masses and widths, radii).
The selection efficiency is a source of systematic uncertainty since the detector simulation may not be perfect.To reduce this effect, a reweighting of the simulated sample is applied to match the kinematical distributions of the data before the BDT selection.Its effect on the result is tested by performing again the reweighting after the BDT selection, and repeating the fit.The difference with respect to the nominal fit is taken as a systematic uncertainty.The dependence of the selection efficiencies with respect to the D 0 transverse momentum has also been studied, without a significant effect being found.
There are a couple of sources of systematic uncertainty related to the background description: the finite size of the D 0 sidebands that are used as a proxy, and the technique to describe the background itself.In the first case, the fit is repeated many times with alternative descriptions of the background obtained by resampling the sidebands using a bootstrapping technique [41].In the second case, an alternative technique is used for the weighting of the simulated sample to match the sidebands.The systematic uncertainties are taken from the differences with respect to the nominal fit.
Alternative parametrisations of the resonance shapes are tested.The Gounaris-Sakurai lineshape, used for the ρ(770) 0 resonance, is compared to the RBW lineshape.Alternative solutions for the KK and ππ K-matrices as reported in Refs.[36,37] are compared with the nominal parametrisation.The fit is repeated many times with the values of the masses and widths are randomly drawn from a Gaussian distribution with mean and width taken as their central value and uncertainty.The values of the Blatt-Weisskopf radii for the D 0 , K * (892) 0 and K * (1680) 0 mesons are varied around their nominal value according to a Gaussian distribution with a width defined by the result of the corresponding amplitude fit.The radii of the other resonances are uniformly varied in a range of ±0.2 GeV −1 around their central values, similarly to Refs.[39,42].
The effect of the fit bias is studied by generating samples and fitting them using the same model.Four sets of pseudoexperiments are performed.The first only uses the nominal signal model to test the stability of the fitter.The second uses the signal and background models to test the effect of the background description.The third uses the signal model with the introduction of a wrong flavour assignment of the D 0 candidate.In data, about 0.5% of the D 0 mesons are associated with the wrong muon [10].The effect on the measurement is tested by generating pseudoexperiments in which 0.5% of the sample is generated as D 0 instead of D 0 .Finally, the fourth uses the signal model where a detection asymmetry between the two kaons is introduced.While some intermediate states are CP eigenstates (e.g.D 0 → φ(1020)ρ(770) 0 ), others are not (e.g.D 0 → K 1 (1270) + K − ) and could be affected by a kaon detection asymmetry.The effect is studied with pseudoexperiments in which a momentum-dependent reconstruction asymmetry of the order of 1% is introduced between the K + and K − mesons.
The model building method chosen in this analysis produces one model, which is the best solution given the amplitudes considered and the criteria used for the selection of the amplitudes.By slightly varying this method, alternative models can be produced with similar fit qualities.Three alternative models are considered to assign a systematic uncertainty.The K * (1410) 0 resonance is used as an alternative to the K * (1680) 0 resonance in the first alternative model.In the second alternative model, five additional amplitudes are added to the nominal model to test the effect of the stopping criteria.Finally, the BaBar collaboration has recently observed the decay of the ρ(1450) 0 meson into two kaons [43].This contribution is tested by adding to the nominal model the component D 0 → ρ(1450) 0 ρ(770) 0 in D-wave, with ρ(1450) 0 → K + K − and ρ(770) 0 → π + π − .The results of these three alternative models are described in Appendix A. Pseudoexperiments are performed by generating a signal sample according to the alternative model and fitting it according to the nominal model.For each amplitude, the largest difference with respect to the nominal fit among the three alternative models is assigned as a systematic uncertainty.
The fraction of signal candidates f s is fixed in the fit.To estimate the impact of its statistical uncertainty to the final result, the fit is repeated many times with values of f s sampled from a Gaussian with mean 0.828 and width of 0.003.No significant effect is found.
The breakdown of all the systematic uncertainties on the model parameters and fit fractions is shown in Appendix B for the CP -averaged and the CP -violating fits.
In addition, some cross-checks are performed.In particular it is checked that the choices made during the selection (the cut on m(K S veto and the treatment of multiple candidates) do not bias the results and that resolution effects are negligible.

CP -averaged results
The fit results are summarised in Table 2, which shows the fit parameters and fit fractions of each component in the model.The fit is performed relative to the fixed component  3. The first uncertainty is statistical and the second is systematic.
A few features of the model are worth noting.The D 0 → φ(1020)ρ(1450) 0 and D 0 → K * (1680) 0 K * (892) 0 components appear only in P -wave without their S-and D-wave counterparts, which are also allowed.The a 1 (1260) + resonance is decaying only to φ(1020)π + , while a contribution of K * (892) 0 K + is reported by the PDG [24].Finally, the ρ − ω superposition seems to be different between the decay modes, as shown in Table 4.
The resulting mass and width of the K 1 (1270) + resonance are 1297 ± 1 MeV/c 2 and Table 4: Parameters of the ρ − ω interference for all relevant amplitudes.The first uncertainty is statistical and the second is systematic.These results are compared to the CLEO legacy-data model of Ref. [8].The main components are present in both models.While the D 0 → φ(1020)ρ(770) 0 components are compatible, the D 0 → K * (892) 0 K * (892) 0 components do not have the same hierarchy between the three angular momentum configurations.The contributions of the D 0 → K 1 (1270) ± K ∓ and D 0 → K 1 (1400) ± K ∓ components are different.One possible explanation is that Ref. [8] did not impose CP conservation in the strong decays of the kaon resonances.The only component of the CLEO legacy-data model that is not found in this analysis is D 0 → K * (1680) + K − , although it is included in the list of potential amplitudes.Instead, other amplitudes with small fit fractions are uncovered in this analysis, as a consequence of the more than 50 times larger dataset.

CP -violation results
For the CP -violation fit, the dataset is split into two subsets according to the charge of the muon to separate the D 0 and CP -transformed D 0 decays.The CP -violation fit described in Sec.4.5 is applied to these two samples.Table 5 shows the resulting CP -violation parameters.The average moduli and phases are not shown in the table as they are identical to the moduli and phases from the CP -averaged fit in Table 2.
All the asymmetry parameters are compatible with zero.The most significant deviation, observed for the phase difference for the D 0 → [φ(1020)ρ(1450) 0 ] L=1 component, corresponds to a 2.8σ statistical fluctuation.To check how likely such a deviation would be in absence of CP violation, the fit is repeated many times, where the data is randomly split instead of splitting it according to the flavour of the D candidate.The largest deviation among all the asymmetry parameters exceeds 2.8σ in 35% of the fits, confirming that the deviation observed in the CP -violation fit is not significant.

Conclusion
An amplitude analysis of the Cabibbo-suppressed decay mode D 0 → K + K − π + π − is performed.The resulting amplitude model provides the most precise description of this decay to date and is used to perform a search for CP violation.
More than 25 decay amplitudes of the D 0 meson have been identified.The most abundant being D 0 → [φ(1020)(ρ − ω) 0 ] L=0 , followed by , and D 0 → K 1 (1270) + K − , all together representing about 80% of the total decay rate (neglecting interference).This model confirms the main findings of Ref. [8] and provides an improved description of the data.In particular, a ρ − ω interference is found that does not allow to treat the two resonances separately and the contribution of the KK, ππ and Kπ S-waves are studied.
For each component of the model, CP asymmetries related to the amplitude modulus, amplitude phase and fit fraction are measured with a total uncertainty ranging from 1% to 15%, dominated by the statistical uncertainty.At this level of sensitivity, no effect of CP violation is found.This is expected from SM predictions [4] and large effects from beyond the SM processes are ruled out.It is interesting to report that most of the systematic uncertainties of the CP -averaged fit are marginally affecting the CP -violation fit, resulting in a much smaller systematic uncertainty.Finally, the systematic uncertainties of the CP -violation measurements are expected to scale down with luminosity, since most of them are estimated with inputs from data.
Table 5: CP -violation parameters fitted simultaneously to the D 0 and (CP -transformed) D 0 samples.The first uncertainty is statistical and the second is systematic.

Appendices A Alternative models
Three alternative models are considered to assign a systematic uncertainty.The results of the corresponding fits are listed in Tables 6-8.
Table 6: Modulus and phase of the fit parameters along with the fit fractions of the amplitudes included in the alternative model using the K * (1410) 0 instead of the K * (1680) 0 resonance.

B Systematic uncertainties
All the systematic uncertainties from all the sources considered on every amplitudes for the nominal and CP -violation fits are listed in Tables 9-12.

Figure 1 :
Figure 1: Mass distribution of the D 0 → K + K − π + π − candidates after the final selection, with fit result superimposed.The top plot shows the normalised residuals.

Figure 2 :
Figure 2: Definition of the helicity angles θ K and θ π , and the decay-plane angle φ.

Table 2 :
Modulus and phase of the fit parameters along with the fit fractions of the amplitudes included in the model.The substructures of the three-body resonances are listed in Table

Figure 3 :
Figure3: Distributions of the five CM variables for the selected D 0 and CP -transformed D 0 candidates (black points with error bars).The results of the five-dimensional amplitude fit is superimposed with the signal model (dashed blue), the background model (dotted green) and the total fit function (plain red).The plot on top of each distribution shows the normalised residuals, where the error is defined as the quadratic sum of the statistical uncertainties of the data and simulated samples.

Table 1 :
Resonances considered in the analysis, classified according to their spin-parity J P and decay products.

Table 9 :
Statistical and systematic uncertainties (in %) on the fit fractions.Values smaller than 0.0005% are displayed as "0.000".The sources of systematic uncertainty are described in the text in the same order as shown in this table.

Table 10 :
Statistical and systematic uncertainties on the fit parameters shown for all floating components.Values smaller than 0.0005 are displayed as "0.000".The sources of systematic uncertainty are described in the text in the same order as shown in this table.For each amplitude, the first value quoted is the modulus and the second is the phase of the complex fit parameter.

Table 11 :
Statistical and systematic uncertainties (in %) on A F k .Values smaller than 0.0005% are displayed as "0.000".The sources of systematic uncertainty are described in the text in the same order as shown in this table.

Table 12 :
Statistical and systematic uncertainties (in %) on the CP -violation parameters shown for all floating components.Values smaller than 0.0005% are displayed as "0.000".The sources of systematic uncertainty are described in the text in the same order as shown in this table.For each amplitude, the first value quoted is the modulus asymmetry and the second is the phase difference.