Search for CP violation in D + ( s ) → K − K + K + decays

A search for direct CP violation in the Cabibbo-suppressed decay D + s → K − K + K + and in the doubly Cabibbo-suppressed decay D + → K − K + K + is reported. The analysis is performed with data collected by the LHCb experiment in proton-proton collisions at a centre-of-mass energy of 13 TeV corresponding to an integrated luminosity of 5.6 fb − 1 . The search is conducted by comparing the D +( s ) and D − ( s ) Dalitz-plot distributions through a model-independent binned technique, based on ﬁts to the K − K + K + invariant-mass distributions, with a total of 0.97 (1.27) million D + s ( D + ) signal candidates. The results are given as p -values for the hypothesis of CP conservation and are found to be 13.3% for the D + s → K − K + K + decay and 31.6% for the D + → K − K + K + decay. No evidence for CP violation is observed in these decays.

To date there is only one observation of CPV in the charm sector, through the difference of CP asymmetries in D 0 → K − K + and D 0 → π − π + decays [19], which is mainly a direct CPV observable [20,21]. New evidence indicates the main source for this CPV effect to be from the D 0 → π − π + channel [22]. Although the measured CP asymmetry, at the 0.1% level, is consistent with some estimates within the SM [11][12][13][14] other estimates indicate that BSM contributions are required to explain this result [23][24][25]. It is thus important to extend the studies to different charm-hadron species decaying into a broad range of final states, including not only CS but also DCS decays [4].
Direct CPV occurs when a given final state is produced through amplitudes with different weak phases while the presence of different strong phases is also required [26]. Three-body charm decays are particularly interesting in this context: final states are reached mainly through resonances, providing an important source of strong-phase differences that vary across the two-dimensional Dalitz plot [27]. This feature can enhance the sensitivity to CP asymmetries in localised regions of the phase space when compared to global charge asymmetries (i.e. integrated over the Dalitz plot) [28][29][30][31] -as observed in charmless B + decays [32].
This paper describes a model-independent search for direct CPV in the Dalitz plots of the CS D + s → K − K + K + and DCS D + → K − K + K + decays 1 with data collected by the LHCb experiment from 2016 to 2018 in proton-proton (pp) collisions at a centre-ofmass energy √ s = 13 TeV corresponding to an integrated luminosity of 5.6 fb −1 . The applied search method is a variation of the original Miranda technique [33,34], which has been used in several studies [35][36][37][38][39]. The method involves dividing the Dalitz plot in two-dimensional bins and computing, for each bin, the significance of the difference in the numbers of D + (s) candidates and D − (s) candidates, where the latter is corrected for global charge asymmetry. A two-sample χ 2 test is then performed on the D + (s) and D − (s) samples. In the original Miranda approach, these numbers are obtained by counting the candidates within a defined K − K + K + invariant-mass signal region, including background. In the approach used here, instead, the yields of D + (s) and D − (s) signal candidates are obtained in each Dalitz-plot bin by a fit to the K − K + K + invariant-mass spectrum of candidates lying in that bin, thus background candidates are not counted. Given the number of degrees of freedom, the probability (p-value) of observing a χ 2 greater than the measured value is computed under the null hypothesis, which in this case is CP conservation. The method is insensitive to global sources of charge asymmetries, but relies on the absence of residual local nuisance asymmetries, for instance from production and detection. This assumption is validated using simulation, by studying the background in the K − K + K + samples, and by analysing large samples of the Cabibbo-favoured modes D + s → K − K + π + and D + → K − π + π + for which CPV is not allowed.

LHCb detector
The LHCb detector [40,41] 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 magnetic field deflects oppositely charged particles in opposite directions which can lead to detection asymmetries. Periodically reversing the magnetic field polarity throughout data-taking cancels most of this effect. The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV. 2 The minimum distance of a track to a primary pp collision 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. 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, and 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, which consists of a hardware stage followed by a two-level software stage. In between the two software stages, an alignment and calibration of the detector is performed in near real-time and their results are used in the trigger [42]. This offers the opportunity to perform physics analyses directly using candidates reconstructed in the trigger [43,44] which the present analysis exploits. Only the reconstructed quantities of triggered candidates are stored, which reduces the event size by an order of magnitude. Trigger decisions at the hardware stage are associated with reconstructed particles, and the decisions with associations are available in the offline selection. Hence requirements are made on the trigger decision on whether it was due to the signal candidates, due to other particles in the pp collision, or a combination of both.
Simulation is used to model the shape of the signal K − K + K + mass distribution accounting for the effects of the detector acceptance, reconstruction and selection criteria. In the simulation, pp collisions are generated using Pythia [45,46] with a specific LHCb configuration [47]. Decays of unstable particles are described by EvtGen [48], in which final-state radiation is generated using Photos [49]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [50] as described in Ref. [51]. The underlying pp interaction is reused multiple times, with an independently generated signal decay for each [52].

Candidate selection and final data samples
The D + → K − K + K + and D + s → K − K + K + decay candidates are selected online by a dedicated software trigger. Three charged particles identified as kaons are combined to form a good-quality decay vertex, detached from any PV. The PV with the smallest value of χ 2 IP is associated to the decay candidate, where χ 2 IP is defined as the difference in the vertex-fit χ 2 of the PV reconstructed with and without the particle under consideration, in this case the D + (s) candidate. Further requirements are applied on the D + (s) decay time; on the angle between the reconstructed D + (s) momentum vector and the vector connecting the PV to the decay vertex; on the χ 2 of the D + (s) decay vertex fit; and on the momentum, the transverse momentum and the χ 2 IP of the D + (s) candidate and its decay products. The invariant masses of the D + and D + s candidates are required to be within the interval 1805-1935 MeV and 1905-2035 MeV, respectively.
In the offline selection, further particle identification (PID) requirements are applied to the decay products to reduce the cross-feed from other charm-hadron decays. For the D + s → K − K + K + channel, the most significant contribution comes from the D + → K − K + π + decay, which is reduced by a stringent PID requirement on the two samecharge kaons. This also reduces contributions from Λ + c → K − K + p and Λ + c → K − π + p decays. For the D + → K − K + K + channel, the only significant cross-feeds are the aforementioned Λ + c decays, which are sufficiently reduced by the PID requirements that maximise the statistical significance for the signal.
A small fraction of D + (s) candidates is rejected when two tracks or segments of tracks are duplicates of each other, resulting in reconstructed same-charge kaon candidates with very similar momenta. Fiducial requirements are applied in order to remove small regions of the phase space where large charge asymmetries are observed, typically caused by low-momentum tracks being swept out of the detector acceptance by the magnetic field.
The remaining background, mostly due to random combinations of tracks, is reduced by using a boosted decision tree (BDT) classifier [53,54], as implemented in the TMVA toolkit [55,56]. Only quantities related to the three-track combinations are used as inputs for the BDT classifier. In addition to quantities used in the online selection, such as flight direction, the quantities IP, χ 2 IP , the distance between the PV and the D + (s) decay vertex and its significance are also used. Approximately 3% of the data is used to train the algorithm and discarded for the rest of the analysis. The signal and the background samples are obtained from this fraction of the data via the sPlot technique [57] using the K − K + K + invariant mass as discriminating variable. Separate BDT classifiers are trained for the D + s → K − K + K + and D + → K − K + K + samples. The requirements on the BDT outputs are chosen to maximise the statistical significance of the signal yields.
The K − K + K + invariant-mass distributions for the final samples are shown in Fig. 1, with the results of maximum-likelihood fits overlaid. The signal shape is described by a Gaussian and two Crystal-Ball [58] functions (CB) with a common mean. Each CB has a width defined as a factor times the Gaussian width. The CB tails are on opposite sides. All the parameters describing the signal shape, apart from the Gaussian width and mean, are fixed to the values obtained from simulation, including the relative fraction of each component. The background is parameterised by a third-order Bernstein polynomial. The selection results in signal yields of 0.97 million D + s → K − K + K + decays and 1.27 million D + → K − K + K + decays. Roughly 1% of the events contain multiple signal candidates. The purity of the samples in a region comprising 95% of the signal candidates is about 64% and 78% for D + s and D + decays, respectively. The Dalitz plots are defined in terms of the variables s high and s low , which represent the higher and lower values of the squared invariant masses formed by the two K − K + combinations. The momenta used to compute these quantities are obtained from a kinematic fit [59] in which the invariant mass of the reconstructed candidates is constrained to the known D + (s) mass [2]. These Dalitz plots can be seen in Fig. 2 for candidates within the K − K + K + invariant-mass region comprising 95% of the signal candidates. For both decay modes the contribution of the φ(1020)K + channel is visible. An amplitude analysis was performed recently for the D + → K − K + K + decay [60] where the f 0 (980)K + and f 0 (1370)K + channels were also found to contribute. To date, no amplitude analysis exists for the D + s → K − K + K + decay, but its Dalitz distribution seems to follow qualitatively the same pattern of that of the D + decay: a clear φ(1020) signature plus a rather smooth distribution elsewhere.

Method
The binned model-independent technique used in this analysis compares the Dalitz-plot distributions for particles and antiparticles, and it is a variation of the original Miranda technique [33,34]. For each Dalitz-plot bin, the local CP observable S CP is defined as 3 yield to that of D + (s) . The method is therefore only sensitive to local CP asymmetries, while global effects such as production and detection asymmetries, if uniform across the Dalitz plot, are corrected for. In the absence of CPV, the values of S i CP follow a standard normal distribution. A two-sample χ 2 test is conducted, with χ 2 = (S i CP ) 2 , and with the number of degrees of freedom (ndof) being the number of bins minus one (due to the constraint of the α normalisation). The resulting p-value from this test is defined as the probability of obtaining a χ 2 that is at least as high as the value observed, under the assumption of CP conservation (null hypothesis). The criterion used is that CP violation is observed for p-values less than 3 × 10 −7 .
The original Miranda method relies on counting the number of particle and antiparticle candidates in each bin, including contributions from both signal and background, with corresponding uncertainties given by their square roots. In this paper, a novel approach is introduced, in which the test is performed by obtaining N i (D ± (s) ) and corresponding uncertainties through fits to the K ∓ K ± K ± invariant-mass distribution of the candidates in each bin. By doing so, the effect of the background contribution to the calculation of S CP is removed. This is particularly relevant for large samples with a significant level of background: when a source of global charge asymmetry (such as that from production effects) affects the signal and the background differently, the α factor calculated using the original Miranda method, which includes contributions from both signal and background, may introduce a bias as shown below.  The binning scheme utilised in this analysis consists of 21 bins for each decay channel, 4 chosen such that the number of signal candidates varies within roughly a factor of two among all bins. The division of the bins around s low ∼ 1 GeV 2 is defined to enhance the sensitivity to CP asymmetries that may change sign when crossing the φ(1020) resonance mass or the node of its angular distribution, for instance due to S-and P-wave interference [61,62]. The binning maps can be seen in Fig. 2. An alternative binning with 50 bins is used for cross-checks. For the K − K + K + invariant-mass fits in each Dalitz-plot bin, the signal and background models are those described in Sect. 3, where the signal shape parameters are fixed to values obtained from simulation in each bin. Fits to the K − K + K + and K + K − K − mass spectra in each bin are performed independently. Figures 3 and 4 show examples of the mass fits for three representative bins.
Studies are performed to validate the fit-per-bin method. Pseudoexperiments are  generated with K − K + K + mass distributions following those obtained in data for D + s and D − s decays, for the 21 Dalitz-plot bins. The total number of events is varied according to a Poisson distribution based on the yields obtained from fitting the samples in each bin without charge separation. The number of candidates, for both signal and background, is scaled by the relative fraction of D + s and D − s candidates in the the full data sample, in order to introduce an integrated asymmetry. With this test, it is possible to compare the response of the original Miranda method with the fit-per-bin method. The distributions of the p-values for the D + s → K − K + K + decay following the two methodologies can be seen in Fig. 5. For the current dataset size and level of background, the original Miranda method introduces a bias towards lower p-values whereas the fit-per-bin method leads to a uniform distribution as expected for the null-hypothesis. The fit-per-bin distribution of p-values is also seen to be uniform for D + → K − K + K + pseudoexperiments.
A cross-check of the fit-per-bin method is also performed using data, where the final samples are divided randomly in two parts according to the overall proportions of D + (s) and D − (s) but with no charge separation. The method is applied for both channels and using the baseline 21-bin scheme and also the alternative 50-bin scheme. The resulting p-values lie between 30% and 91%, compatible with the null-hypothesis. The approximate sensitivity of the current data samples to CP -violating effects is studied using simulated pseudoexperiments with similar resonant structure [60] and number of signal and background candidates as in data. In these pseudoexperiments, phases or relative magnitudes are set to different values in the key resonant states, such as φ(1020)K + and f 0 (980)K + , between the D + (s) and D − (s) amplitudes. It is found that the p-value threshold of 3 × 10 −7 can be achieved for a difference of 3 • to 7 • in phase or 3% to 7% in relative magnitude depending on the resonant state for which these effects are introduced in D + (s) and D − (s) decays.

Sensitivity to nuisance asymmetries
The sensitivity to possible nuisance localised charge asymmetries induced by detection efficiency or production rate is assessed through studies using fully simulated samples, the K − K + K + background candidates in data, and Cabibbo-favoured control samples.
Using simulation samples, about 2.7 and 5.4 million D + s → K − K + K + and D + → K − K + K + candidates, respectively, are retained after applying the same selection requirements as those in data. Since there is no background and the samples were generated with a flat distribution in phase space, the original Miranda method is used with a 21 samesize binning scheme. The p-values found are 37.7% and 76.4% for the D + s → K − K + K + and D + → K − K + K + decays, respectively, showing no induced asymmetries from the detector responses, according to simulation.
For the background study, there are in total about 3.4 and 2.5 million candidates in the D + s and D + samples, respectively. The consistency with no charge asymmetry is tested using the fit-per-bin method described in Sect. 4, where the signal yields and uncertainties are replaced with those of the background from the fits per bin. The resulting p-values are 46.0% and 75.7% for the D + s and D + backgrounds, respectively, which indicate no effect of local charge symmetries.
The Cabibbo-favoured D + s → K − K + π + and D + → K − π + π + decays are used as control channels. These modes have much larger data samples than the analysed signal samples. To test whether local nuisance asymmetries are present in the data, the control samples are divided into multiple subsets of approximately five times the size of the signal channels. There are 45 samples of D + s → K − K + π + decays and 85 of D + → K − π + π + decays with over 90% purity. Since the selection criteria are different for the D + s → K − K + K + and D + → K − K + K + samples, the control-channel samples are studied applying both sets of criteria (except the PID requirements for the pion candidates). To perform the search for local asymmetries, four binning schemes are tested: two uniform grids, of 5 × 5 and 8 × 8, and two adaptive binnings such that all bins have the same population, with 21 and 50 bins in total. Due to the high purity of the samples, the original Miranda technique is used, i.e. simple counting to obtain the number of D + s and D + candidates, respectively, in the m(K − K + π + ) region of  MeV and in the m(K − π + π + ) region of 1850-1890 MeV. The resulting p-value distributions for both control samples under the different selection criteria and binning schemes are found to be uniformly distributed, with no value incompatible with the null hypothesis. These results are thus consistent with the absence of local nuisance charge asymmetries. Figures 6 and 7 show examples of the S CP distribution for one control-mode subset and the overall p-value distribution for the adaptive 21-bin scheme for the D + s → K − K + π + and the D + → K − π + π + channels where the selection criteria of the D + s → K − K + K + channel are applied to both. Values of S CP distribution for one subset of D + → K − π + π + data with D + s → K − K + K + selection criteria. Right: Overall p-value distribution for the 85 D + → K − π + π + subsets using the scheme of 21 adaptive bins.

Results
Searches for CPV in the D + (s) → K − K + K + decays are performed with the fit-per-bin strategy described in Sect. 4. The total signal yields for each charge and the values of α for the D + s → K − K + K + and D + → K − K + K + channels are displayed in Table 1. The distributions of S CP for both decays are shown in Fig. 8. The resulting p-values are 13.3% for the decay D + s → K − K + K + and 31.6% for the decay D + → K − K + K + . The results are consistent with the hypothesis of no localised CP violation in either channel.
The dependence of the results on the signal and background models is studied. The signal shape is changed to two Gaussians plus two CB functions, the background is changed to a third-order Chebyshev polynomial, and combinations of these changes are also tried. Additionally, instead of independent fits, simultaneous fits where the background shape parameters are shared among D + (s) and D − (s) samples are also tested. Under these variations, the p-values for the D + s → K − K + K + decay range from 1.7% to 20.2%, while for D + → K − K + K + they range from 3.2% to 48.5%.
An alternative binning scheme is also tested, with 50 bins, using the baseline models for signal and background. A finer binning may allow a higher sensitivity in localised regions, though it could lead to an overall loss of statistical sensitivity. The S CP distributions for this alternative binning are shown in Fig. 9. The resulting p-values are 10.1% and 8.2%    for the D + s → K − K + K + and D + → K − K + K + channels, respectively. Altogether the results consistently indicate no evidence for CP violation in the current samples for D + s → K − K + K + and D + → K − K + K + decays.

Conclusions
This paper presents a search for CP violation in the Cabibbo-suppressed decay D + s → K − K + K + and in the doubly Cabibbo-suppressed decay D + → K − K + K + us-ing pp collison data at √ s = 13 TeV and corresponding to an integrated luminosity of 5.6 fb −1 . The study is performed with samples containing about 0.97 (1.27) million D + s (D + ) → K − K + K + decays. A novel technique is used: the fit-per-bin method, which is a variation of the original Miranda technique, where the value of S CP is calculated in each Dalitz plot bin using the D + (s) and D − (s) signal yields estimated from K − K + K + mass-invariant fits. The fit-per-bin method is found to be robust for cases with very large signal samples and significant background.
The results are given as p-values with respect to the null-hypothesis of CP conservation and are found to be 13.3% for the D + s → K − K + K + channel and 31.6% for the D + → K − K + K + channel. No evidence for CP violation is found. Cross-checks by varying the invariant-mass fit models and binning schemes give consistent results. This is the first search for CP violation in the Cabibbo-suppressed channel D + s → K − K + K + and in the doubly Cabibbo-suppressed channel D + → K − K + K + . The work presented here contributes to a wider effort of comprehensive studies of CP violation in charm decays.