Searches for 25 rare and forbidden decays of $D^+$ and $D_s^+$ mesons

A search is performed for rare and forbidden charm decays of the form $D_{(s)}^+ \to h^\pm \ell^+ \ell^{(\prime)\mp}$, where $h^\pm$ is a pion or kaon and $\ell^{(')\pm}$ is an electron or muon. The measurements are performed using proton-proton collision data, corresponding to an integrated luminosity of $1.6\text{fb}^{-1}$, collected by the LHCb experiment in 2016. No evidence is observed for the 25 decay modes that are investigated and $90\%$ confidence level limits on the branching fractions are set between $1.4\times10^{-8}$ and $6.4\times10^{-6}$. In most cases, these results represent an improvement on existing limits by one to two orders of magnitude.


Introduction
The search for rare and forbidden decays in flavour physics constitutes an important test of the Standard Model (SM) and provides a window to new physics. In the charm sector, decays of the form D + → h ± + ( )∓ and D + s → h ± + ( )∓ , where h ± is a charged pion or kaon and l ( )± is an electron or muon, are among such processes. 1 Searches are reported for 25 decay modes, with common experimental approaches applied to all channels.
The physics processes of these channels vary. Four of the decay channels (D + → π + e + e − , D + → π + µ + µ − , D + s → K + e + e − , D + s → K + µ + µ − ) involve flavour-changing neutral current (FCNC) transitions. These processes are rare within the Standard Model as only weak annihilation processes can occur at tree level. At the loop level, FCNC transitions are suppressed by the Glashow-Iliopoulos--Maiani (GIM) mechanism [1] but are well established in both B-meson and K-meson decays [2]. FCNC B-meson decays have received significant attention in recent years: the process B 0 s → µ + µ − was observed [3,4]; and experimental observations in analyses of b → s + − transitions are in tension with the SM predictions (see Ref. [5] and references therein). The GIM cancellation in the D-meson system is stronger than in the B-meson system, leading to short-distance SM branching fractions of O(10 −12 ) [6]. Significant long-distance contributions occur from dilepton resonances and these are expected to dominate the SM contribution across the full dilepton invariant-mass-squared distribution [6,7].
The potential of c → u + − processes to constrain new physics has been discussed in the literature [6][7][8][9]. Model-independent constraints on Wilson coefficients in effective field theory have been considered, as have specific classes of models. Recently particular attention has been paid to leptoquark models, due to their potential relevance to the anomalies on b → s + − processes [10][11][12][13][14].
The analysis reported here is performed on a data sample corresponding to an integrated luminosity of 1.6 fb −1 collected by the LHCb experiment in 2016. Four resonant decays are used for calibration and normalisation. Non-resonant decays with a pion and two same-sign or opposite-sign muons in the final state have previously been searched for by LHCb with no evidence observed [15]. The best limits in the other channels are measured by the BaBar collaboration [16], which has studied all channels, or by the CLEO collaboration, from searches in decay modes with dielectrons in the final state [17].

Triggering, reconstruction and selection
The online event selection is performed by a trigger, which consists of a hardware stage, based on information from the calorimeter and muon systems, 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 the results are used in the trigger. The same alignment and calibration information is propagated to the offline reconstruction, ensuring consistent and high-quality particle identification (PID) performance between the trigger and offline software.
A hardware trigger selection is made for all decay channels based on observing a deposit with high transverse energy in the hadronic calorimeter. Depending on the leptons in the final state, this is supplemented by a selection of events with high transverse-energy clusters in the electromagnetic calorimeter and high transverse-momentum muons. At the first software level, inclusive multivariate selections are used to select charged tracks which do not originate from any PV.
The second software level builds D + (s) candidates using exclusive selections for each combination of final-state particles. Only events containing a single reconstructed PV are used. All three tracks are required to be incompatible with coming from the PV, and to have been reconstructed with good fit quality, have a transverse momentum exceeding 300 MeV and a momentum exceeding 2000 MeV. The tracks are required to form a D + (s) secondary vertex, with good fit quality and all pairs of tracks intercepting within 150 µm. The angle between the reconstructed D + (s) and the vector connecting the PV to the decay vertex of the D + (s) candidate (direction angle) must be within 14 mrad. Loose PID requirements are made on all final-state tracks. Eight of the signal channels in this analysis can proceed via decays through intermediate resonances (η, ρ 0 , ω and φ). These are removed by vetoing the region [525 MeV, 1250 MeV] in dilepton mass. The is further selected by requiring the invariant mass of the dilepton pair be within ±20 ( +40 −100 ) MeV of the φ mass [2]. A bremsstrahlung reconstruction procedure is used to correct the momentum of electron candidates.
A classifier trained using XGBoost [26], is used to further distinguish between signal and background originating from random combinations of tracks. A single classifier is trained for each final state, with signal being represented by an equally weighted mixture of simulated D + and D + s events. This approach is found to give equivalent performance when compared to using a separate classifier for each signal meson. Background candidates for the training are provided from a data sample where all three final-state particles have the same electric charge. The classifier is trained using the fit quality of the primary and secondary vertex, the D + (s) pseudorapidity, the D + (s) flight distance, the difference in χ 2 of the PV reconstructed with and without each final-state track, the momentum of each final-state track, the maximum distance of closest approach between all final-state tracks, the direction angle and the reconstructed proper lifetime of the D + (s) . To guarantee that no other tracks can be associated to the secondary vertex, an isolation variable, A p T , is used that considers the imbalance of p T of nearby tracks compared to that of the D + (s) candidate, where p T (D + (s) ) is the p T of the D + (s) meson and ( − → p ) T is the transverse component of the vector sum of all charged particles momenta within a cone around the candidate, excluding the three signal tracks. The cone is defined by a circle of radius 2.0 in the plane of pseudorapidity and azimuthal angle, measured in radians, around the D + (s) candidate direction. The signal D + (s) candidates tend to be more isolated than the combinatorial background and show on average greater values of A p T . The final selection was optimised using a grid search in the thresholds applied to the classifier output and the PID variables.

Invariant-mass distributions
The signal yields are determined using maximum-likelihood fits to the invariant-mass distributions in each of the D + (s) → h ± + ( )∓ final states, including the resonant channels. The fit is performed in the mass range of 1802 to 2050 MeV for all final states where both the D + and D + s decay channels are analysed. In the final states where only the D + s contribution is analysed, the fit is performed in the range of 1926 to 2050 MeV. The invariant-mass distributions are shown in Figs. 1-5, with fits indicated. No evidence is observed in any of the 25 signal channels. The yields obtained from the fit to the normalisation channels can be found in Table 1.
The mass distributions for the signal contributions are obtained by modelling the distributions in simulated events with a kernel density estimation (KDE) technique [27].   The same technique is applied to the D + (s) → (φ → e − e + ) π + channels. In the D + (s) → (φ → µ − µ + ) π + channels, where a clean separation from background can be made, a fit with a double Gaussian distribution is performed to the data, with the two Gaussian distributions sharing a common mean.
In decay channels with electrons, the observed mass distribution is dependent on the number of reconstructed bremsstrahlung photons. For final states with one electron, separate shapes are extracted from simulation for candidates with no bremsstrahlung and where one or more reconstructed photon candidates are associated with the electron. For final states with two electrons, separate shapes are used for candidates with no bremsstrahlung, one photon per D + (s) candidate, and two or more photons per D + (s)

candidate.
Background arises from random combinations which vary smoothly across the fit range, and from misidentified and partially reconstructed backgrounds that may peak in the fitted mass range. The combinatorial backgrounds are fitted with an exponential distribution where the analysed signal channel contains dimuons, while third-order Chebyshev polynomials are used in the cases where the signal channel contains one or two electrons. The form used for modelling these combinatorial backgrounds is determined using a data sample where all three final-state tracks have the same charge.
Background from decays involving leptons and neutrinos are found to be negligible for this analysis, with only three-body D + and D + s hadronic decays contributing to the signal samples. The backgrounds that affect a specific signal channel depend on whether the signal-channel hadron is a pion or a kaon and whether the leptons are of the same or opposite charge. Simulation samples, which do not use a description of the interaction of the particles with the detector material, are generated [28], and are fitted using the same KDE technique as the signal.
The modelling of the backgrounds is validated using data samples where the nominal kinematic selection is applied, but with the particle identification requirements reversed and their values modified from those in the signal selection to enhance the background.
The distributions in these enhanced background samples are then fitted with the shapes obtained from the KDE fit to the simplified simulation plus a combinatorial-background description. The four final states of the form D + (s) → K + + ( )− have significant contributions from D + (s) → K + π + π − decays, and the seven final states of the form D + (s) → π + ± ( )∓ have significant contributions from D + (s) → π + π + π − decays. In these cases, the KDE description is smeared by a Gaussian function with mean and width included as free parameters of the fit. After applying the additional smearing a good description is obtained in the background-enhanced samples.

Branching fraction determination
The signal branching fractions are obtained from fits to the three-body invariant-mass distributions, shown in Figs. 2-5. The branching fractions are defined by where N is the fitted yield in the channel indicated, B is the branching fraction and is the efficiency of the selection criteria. The branching fraction used for D + (s) → (φ → − + ) π + can be found in Table 1 and is obtained by combining the mea- [2]. The systematic uncertainties are described in Section 6 and taken into account with a lognormal distribution in the branching fraction fit. The potential overlap between the D + and D + s signal peaks in the same final state is accounted for by floating the yield of the other meson in each fit and treating this as an additional nuisance parameter when computing the significance of the signal peak.
The effect of discrepancies between data and simulation are reduced by using a reweighting technique that utilises a multivariate classifier. This classifier is used to generate per-event weights using the method outlined in Ref. [29]. The procedure is applied separately for each of the calibration channels and each classifier is trained to distinguish between real and simulated calibration-channel events. Background subtraction of the data is applied using the sPlot technique [30] with a fit to the invariant mass of the calibration channel. The weights obtained from these classifiers are applied to the simulation and then used throughout this analysis. Decay channels with two electrons in the final state use weights from the D + (s) → (φ → e − e + ) π + classifiers and all other decays use the D + (s) → (φ → µ − µ + ) π + classifiers. In all cases the same parent particle (D + (s) ) decay as the signal is used. Candidates/(4 MeV)

LHCb
Non-peaking   Equation (2) includes the ratio of efficiencies between the normalisation channel and the signal channel so that the effects of several systematic uncertainties cancel. The reconstruction and selection efficiencies are obtained from simulation with the aforementioned corrections applied. As the resonant structure in the signal channels is inherently unknown, the efficiency of each signal decay is obtained under the assumption that the decay particles are uniformly distributed across the phase space with the signal contribution extrapolated into any vetoed kinematic regions.
The efficiency with which events pass the PID requirements is obtained using a set of calibration channels. The PID response is sampled [31] as a function of the particle kinematics and event multiplicity in simulation, as described in Ref. [32]. A small correction is made for differences in the track reconstruction efficiency between data and simulation [33].
For channels containing electrons, an additional correction is applied to the ratio of electron to muon reconstruction efficiencies. This is computed under the assumption that the correction for each electron is independent of the other, and thus the square root of the efficiency-corrected signal-yield ratio between the calibration channels D + (s) → (φ → e − e + ) π + and D + (s) → (φ → µ − µ + ) π + in data is taken. The single event sensitivities, i.e. the branching fractions corresponding to a single observed signal event, vary from 3 × 10 −10 to 1 × 10 −8 depending on the channel. To ensure the efficiencies are sufficiently well understood, cross-check measurements are made of the branching fractions of D + (s) → (φ → µ − µ + ) π + and D + (s) → (φ → e − e + ) π + decays. The results obtained before and after applying the offline selection criteria are compared and found to be in agreement.

Systematic uncertainties
The efficiency ratio in Eq. 2 is affected by several sources of systematic uncertainty. The finite size of the simulated signal event samples introduces a systematic uncertainty on the branching fractions varying from 1.4 % to 6.4 %. A summary of the systematic uncertainties assigned for each effect is given in Table 2.
The track-reconstruction efficiency from simulation is corrected using a tag-and-probe technique in data where J/ψ → µ + µ − decays are selected by making requirements on only one muon. This is found to have negligible effect on the efficiency ratio. A 0.8 % systematic uncertainty is assigned for each particle species in the signal decay that differs from that of the D + (s) → (φ → µ − µ + ) π + reference. An additional 1.5 % systematic uncertainty is assigned to channels containing a kaon to account for possible hadronic interaction effects with the detector material [33]. Alternative parametrisations [32] are considered for the sampling of the PID response, but their associated systematic effects have negligible impact compared to the statistical uncertainties of the analysis and other sources of systematic uncertainty. For final states with both a muon and an electron, a further systematic uncertainty of 7.6 % is assigned for the choice of reweighting classifier. This is obtained from the RMS of the change in efficiency of all mixed-lepton final states, using the alternative reweighting classifier.
The systematic uncertainties affecting the signal yield are estimated using an alternative fit model for the normalisation channels. For signal states with two muons, a KDE parametrisation obtained from simulation is used as an alternative for the nominal model describing the D + (s) → (φ → µ − µ + ) π + signal. In all other cases, the differences are most likely to arise from the treatment of bremsstrahlung radiation in simulation. Therefore, an alternative model is generated by fixing the relative yield between the three bremsstrahlung categories of D + (s) → (φ → e − e + ) π + from simulation. The signal yield from the alternative models is compared with the nominal model and the difference is assigned as a systematic uncertainty.
The uncertainty on the yield of D + (s) → (φ → µ − µ + ) π + is dominated by the modelling of the D + (s) → π + π + π − backgrounds in the lower tails of the larger signal components. To account for this, the data set is refitted neglecting these background components, and the change in signal yield is assigned as a systematic uncertainty. The same procedure is applied for D + (s) → (φ → e − e + ) π + with the change in signal yield used as the uncertainty on the electron efficiency correction factor.

Results and conclusions
No significant deviation from the background only hypothesis is found for any of the channels and all limits are within ±2σ of the expected limit. The compatibility of the observed mass-distributions with a signal-plus-background or a background-only hypothesis is evaluated using the CL s method with systematic uncertainties included as sample apply to all channels. The electron efficiency correction, described at the end of Section 5, is denoted by electron .
All values are given in percent as a fractional uncertainty on the signal yield. Channel Simulated signal Track-reconstruction   [34,35]. Upper limits on the branching fractions are determined using the observed distribution of CL s as a function of the assumed branching fraction. The upper limits at 90 % and 95 % confidence level (CL) are given in Table 3 and are shown in Fig. 6.
In conclusion, searches have been made for 25 previously unobserved three-body decays of D + and D + s mesons using 1.6 fb −1 of pp collision data collected by the LHCb experiment during 2016. The decays are of the form D + (s) → h ± + ( )∓ where h is a kaon or pion and ( ) is an electron or muon. No significant deviations from the background-only hypotheses are seen. The 90 % CL limits on the branching fractions vary from 1.4 × 10 −8 to 6.4 × 10 −6 and represent the world's best limits for 23 of these decays. In the majority of the channels the improvement on previous limits is more than an order of magnitude, with the largest improvement being a factor of five hundred in D + s → K − µ + µ + decays.