Search for the $C\!P$-violating strong decays $\eta \to \pi^+\pi^-$ and $\eta^\prime(958) \to \pi^+\pi^-$

A search for the $C\!P$-violating strong decays $\eta \to \pi^+\pi^-$ and $\eta^\prime(958) \to \pi^+\pi^-$ has been performed using approximately $2.5 \times 10^{7}$ events of each of the decays $D^+ \to \pi^+\pi^+\pi^-$ and $D_s^+ \to \pi^+\pi^+\pi^-$, recorded by the LHCb experiment. The data set corresponds to an integrated luminosity of 3.0 fb$^{-1}$ of $pp$ collision data recorded during LHC Run 1 and 0.3 fb$^{-1}$ recorded in Run 2. No evidence is seen for $D^+_{(s)} \to \pi^+ \eta^{(\prime)}$ with $\eta^{(\prime)} \to \pi^+\pi^-$, and upper limits at 90% confidence level are set on the branching fractions, $\mathcal{B}(\eta \to \pi^+\pi^-)<1.6 \times 10^{-5}$ and $\mathcal{B}(\eta^\prime \to \pi^+\pi^-)<1.8 \times 10^{-5}$. The limit for the $\eta$ decay is comparable with the existing one, while that for the $\eta^\prime$ is a factor of three smaller than the previous limit.


Introduction
The strength of CP violation in weak interactions in the quark sector is well below what would be required to serve as an explanation for the observed imbalance between the amounts of matter and antimatter in the universe. The QCD Lagrangian could contain a term, the θ term [1], that would give rise to CP violation in strong interactions; however, no strong CP violation has been observed. The experimental upper limit on the neutron electric dipole moment (nEDM) implies a limit θ 10 −10 [2]. The closeness of the value of θ to zero is seen as a fine-tuning problem, the so-called "strong CP problem". Solutions to the strong CP problem may involve axions [3], extra space-time dimensions [4], massless up quarks [5], string theory [6] or quantum gravity [7].

Detector and simulation
The LHCb detector [11,12] 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, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary pp interaction vertex (PV), the impact parameter, is measured with a resolution of (15 + 29/p T ) µm, where p T is the component of the momentum transverse to the beam, in GeV/c. Different types of charged hadrons are distinguished using information from two ring-imaging Cherenkov detectors. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic calorimeter and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.
The online event selection is performed by a trigger [13], which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction. At the hardware trigger stage, events are required to have a muon with high p T or a hadron, photon or electron with high transverse energy in the calorimeters.
A new scheme for the LHCb software trigger was introduced for LHC Run 2. Alignment and calibration are performed in near real-time [14] and updated constants are made available for the trigger. The same alignment and calibration information is propagated to the offline reconstruction, ensuring high-quality particle identification (PID) and consistent information between the trigger and offline software. The larger timing budget available in the trigger compared to that available in Run 1 also results in the convergence of the online and offline track reconstruction, such that offline performance is achieved in the trigger. The identical performance of the online and offline reconstruction offers the opportunity to perform physics analyses directly using candidates reconstructed in the trigger [15].
In the simulation, pp collisions are generated using Pythia [16] with a specific LHCb configuration [17]. Decays of hadronic particles are described by EvtGen [18], in which final-state radiation is generated using Photos [19]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [20] as described in Ref. [21].

Data samples and outline of analysis method
In the analysis, the decays D + → π + π + π − and D + s → π + π + π − are used to look for the presence of η and η resonances in the π + π − mass spectra, which could come from the known decays D + (s) → π + η ( ) (inclusion of charge-conjugate modes is implied throughout). The data samples comprise about 25 million each of D + → π + π + π − and D + s → π + π + π − decays, from integrated luminosities of 3.0 fb −1 of pp collision data recorded by LHCb in LHC Run 1 and 0.3 fb −1 recorded in 2015 during Run 2. For N (η ( ) ) observed η ( ) signal decays in the π + π − mass spectrum from a total of N (D + (s) ) mesons reconstructed in the π + π + π − final state, the measured branching fraction would be where (η ( ) ) accounts for any variation of efficiency with π + π − mass, as discussed in Sec. 5.2. The values of N (D + (s) ) and N (η ( ) ) and their uncertainties are obtained from fits to the π + π − π + and π + π − mass spectra of the selected D + (s) → π + π + π − candidates; the branching fractions B(D + (s) → π + π + π − ) and B(D + (s) → π + η ( ) ) and their uncertainties are taken from Ref. [22]; and the relative efficiency factors, , are obtained from simulations. Since the analysis starts from a given number of selected D + (s) → π + π + π − decays, there are no normalisation channels. All selections are finalised and expected sensitivities are evaluated before the η and η signal regions in the π + π − mass spectra are examined.

Event selection
The event selection comprises an initial stage in which relatively loose criteria are applied to select samples of candidate D + (s) → π + π + π − decays. A boosted decision tree (BDT) [23] is then used to further suppress backgrounds. Candidate D + (s) → π + π + π − decays are required to have three good quality tracks, each with p T greater than 250 MeV/c, consistent with coming from a vertex that is displaced from any PV in the event. Loose particle identification criteria are applied, requiring the tracks to be consistent with the pion hypothesis. The three-track system is required to have total charge ±e, its invariant mass must be in the range 1820-2020 MeV/c 2 , and its combined momentum vector must be consistent with the direction from a PV to the decay vertex. The invariant mass of opposite-sign candidate pion pairs is required to be in the range 300-1650 MeV/c 2 ; this removes backgrounds where a random pion is associated with a vertex from either a γ → e + e − conversion, in which both electrons are misidentified as pions, or from a D 0 → K − π + decay, where the kaon is misidentified as a pion.
The BDT has six input variables for each of the tracks, together with three variables related to the quality of the decay vertex and the association of the D + (s) candidate with the PV. The track variables are related to track fit quality, particle identification probabilities and the quality of the track association to the decay vertex. The BDT is trained using a sample of 820 000 simulated D + → π + π + π − events for the signal, generated uniformly in phase space, and about 10 7 background candidates obtained from sidebands of width 20 MeV/c 2 on each side of the D + → π + π + π − mass peak in the data.
The selection criteria for the BDT output value and π + π + π − signal mass windows are simultaneously optimised to maximise the statistical significance of the D + (s) signals, N sig / N sig + N bkg , where N sig is the number of D + (s) signal decays and N bkg is the number of background events within the signal mass windows. The BDT selection gives signal efficiencies of 90% while rejecting about 60% of the backgrounds. The optimum mass selection ranges are ±20 MeV/c 2 for both the D + and D + s peaks in Run 1 and ±21 MeV/c 2 for both peaks in Run 2. Figure 1 shows the π + π + π − mass spectra for Runs 1 and 2, after the BDT selection. The discontinuity in the Run 2 spectrum comes from the fact that the trigger has two separate output streams and there are different BDT cuts for D + and D + s . The yield per fb −1 is larger in Run 2 than in Run 1 by a factor 3.3, arising from the larger crosssection [24], and from a higher trigger efficiency for charm. The curves in Fig. 1 show the results of fits to the spectra in which each peak is parametrised by the sum of a double-sided Crystal Ball function [25] and a Gaussian function, while a fourth-order polynomial is used for the combinatorial background. All shape and yield parameters are allowed to vary in the fits. The fits also include components for contributions from D + s → K + π + π − decays, where the kaon is misidentified as a pion, and from D + s → π + π + π − π 0 and D + (s) → π + η ( ) with η ( ) → π + π − γ. The yields for these last components, the shapes for which are obtained from simulation, are found to be small. The total D + (s) → π + π + π − signal yields in the optimised mass windows, summed over Run 1 and Run 2 data, are 2.49 × 10 7 for D + and 2.37 × 10 7 for D + s , with backgrounds of 1.38 × 10 7 and 1.08 × 10 7 , respectively, within the same mass windows. Uncertainties of ±2%, corresponding to the maximum values of the fit residuals, are assigned to each total yield to account for imperfections in the fits to the mass spectra. To improve the π + π − mass resolution, a kinematic fit [26] is performed on the selected D + (s) candidates, with the three tracks constrained to a common vertex, the π + π + π − mass constrained to the known D + (s) mass, and the D + (s) candidate constrained to come from the PV. 5 Limits on the η ( ) → π + π − branching fractions 5.1 Mass spectra for π + π − For each of the η and η resonances there are four separate π + π − mass spectra, from the D + and the D + s for each of Runs 1 and 2. Figures 2 and 3 show the sums of the four π + π − mass spectra for the η and η mass fitting ranges, which are chosen to avoid the peaks from the   of the decays η ( ) → π + π − γ, using the matrix element given in Ref. [27], show that the contributions from these channels are small and do not peak in the fitting ranges. They are therefore considered as part of the background, which is parametrised by a polynomial function (see Sect. 5.3).
Expected signal π + π − mass line shapes for η → π + π − and η → π + π − are obtained from simulations. In both cases a double Gaussian shape is found to describe the signal well, with mass resolutions of 2.3 MeV/c 2 for the η mass region and 3.2 MeV/c 2 for the η region. These results are calibrated by comparing the η mass resolution from the simulation with that for reconstructed K 0 S → π + π − decays from background D + (s) → K 0 S π + events in the data, before the kinematic fits to the D + (s) candidates. The differences, which are 5% in Run 1 and 10% in Run 2, are taken as the systematic uncertainties on the π + π − mass resolution for both the η and η mass ranges.

Relative efficiency as a function of π + π − mass
The relative efficiency factors in Eq. 1 are obtained from simulation. Fully simulated π + π − mass spectra from D + → π + π + π − decays for Run 1 are divided by the generated spectra to give the relative efficiency as a function of the π + π − mass. The efficiency is highest at large π + π − masses, mainly due to the effects of the hardware and software triggers. The relative efficiencies in Run 1 data are found to be (η) = 0.85 ± 0.01 and (η ) = 1.01 ± 0.01, where the uncertainties come from the simulation sample size. The relative efficiencies for Run 2 are found to be statistically compatible with those for Run 1, through a comparison of the π + π − mass spectra from the D + and D + s signal candidates in the data. An additional systematic uncertainty of 2% is assigned to the Run 2 relative efficiencies, corresponding to the maximum difference between the mass spectra.

Sensitivity studies
In order to measure the sensitivity of the analysis, each π + π − mass spectrum is fitted with a fourth-order polynomial, initially with the signal regions excluded. The signal regions are then populated with pseudo data, generated according to the fitted polynomial functions, with Gaussian fluctuations. Each spectrum is then fitted again with the sum of a fourth-order polynomial plus the η ( ) signal function, and Eq. 1 is then used to obtain branching fractions measured with the pseudo data. As expected, these branching fractions are consistent with zero. Expected upper limits on the branching fractions are obtained using the CL s method [28]. In each case, CL s values are obtained using the products of the likelihood functions for the four individual spectra. Systematic uncertainties are included, but have no effect on the results, which are shown in Fig. 4 for the η and in Fig. 5 for the η . Expected limits at 90% CL are B (η → π + π − ) < 2.0 × 10 −5 and B (η → π + π − ) < 1.8 × 10 −5 .

Observed limits on the branching fractions
The procedures outlined above are then applied to the observed mass spectra, i.e. with the pseudo data in the signal ranges replaced by the observed data. The sums of the fits to the four spectra for the η and η are shown as the solid curves in Figs. 2 and 3. The results are consistent among the four mass spectra for each meson. Weighted average branching fractions are measured to be B (η → π + π − ) = (−1.1 ± 1.8) × 10 −5 and B (η → π + π − ) = (0.8 ± 1.6) × 10 −5 . Although the simple, unweighted sum of the fits to the π + π − mass spectra in Fig. 2 shows a small, but insignificant, positive yield for the η, the weighted average branching fraction B (η → π + π − ) is dominated by a negative value in the Run 1 D + s sample.  Since there is no evidence for any signal, the CL s method is used, as for the pseudo data, to obtain observed upper limits on the branching fractions. Figures 4 and 5 show the observed CL s values as functions of the branching fractions. Limits obtained are B η → π + π − < 1.6 × 10 −5 , B η → π + π − < 1.8 × 10 −5 , both at 90% CL, in good agreement with the expected limits.

Conclusions
A new method is introduced to search for the decays η → π + π − and η (958) → π + π − , which would violate CP symmetry in the strong interaction. The method relies on the copious production of charm mesons at LHCb, and will improve in sensitivity as more data are collected at the LHC. With the LHC Run 1 data and data from the first year of Run 2, the limit obtained on the branching fraction for the decay η → π + π − is comparable to the existing limit, while that for η → π + π − is a factor three better than the previous limit.