Searches for 25 rare and forbidden decays of D+ and Ds+\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {D}_s^{+} $$\end{document} mesons

A search is performed for rare and forbidden charm decays of the form Ds+→h±ℓ+ℓ′∓\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {D}_{(s)}^{+}\to {h}^{\pm }{\mathrm{\ell}}^{+}{\mathrm{\ell}}^{\left(\prime \right)\mp } $$\end{document}, where h± is a pion or kaon and ℓ(′)± is an electron or muon. The measurements are performed using proton-proton collision data, corresponding to an integrated luminosity of 1.6 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 × 10−8 and 6.4 × 10−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 potential of c → u + − processes to constrain new physics has been discussed in the literature [8][9][10][11]. 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 [12][13][14][15][16].
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 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 [17]. The best limits in the other channels are measured by the BaBar collaboration [18], which has studied all channels, or by the CLEO collaboration, from searches in decay modes with dielectrons in the final state [19].

Detector and simulation
The LHCb detector [20,21] 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 silicon-strip 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 the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV. The minimum distance in the plane transverse to the beam 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, an electromagnetic and a hadronic calorimeter. Muons are identified by a system composed of alternating layers of iron and multiwire proportional chambers.

JHEP06(2021)044
Simulation is required to model the effects of the detector acceptance and the selection requirements. In the simulation, pp collisions are generated using Pythia [22,23] with a specific LHCb configuration [24]. Decays of unstable particles are described by EvtGen [25], in which final-state radiation is generated using Photos [26]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [27,28] as described in ref. [29].

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 is further selected by requiring the invariant mass of the dimuon (dielectron) 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 [30] 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 an identically selected data sample except all three final--3 -JHEP06(2021)044 state particles are required to 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 three-dimensional 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, 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 figures 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 [31]. The same technique is applied to the D + (s) → φ → e − e + π + channels. In the 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 function used for modelling these combinatorial backgrounds is determined using a data sample where all three final-state tracks have the same charge and any parameters are left floating in the final fit.

Non-peaking
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 [32], 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. Candidates/(4 MeV)

Branching fraction determination
The signal branching fractions are obtained from fits to the three-body invariant-mass distributions, shown in figures 2-5. The branching fractions are defined by

JHEP06(2021)044
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 measurements of D + (s) → φ → K + K − π + , φ → K + K − and φ → + − from ref. [2]. The systematic uncertainties are described in section 6 and taken into account with a log-normal 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.
Data and simulation are found to be within good agreeement and the quality of the modelling for the three-body invariant mass can be seen in figure 1. The effect of discrepancies between data and simulation are further 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. [33]. 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 [34] 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. Equation (5.1) 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 [35] as a function of the particle kinematics and event multiplicity in simulation, as described in ref. [36]. A small correction is made for differences in the track reconstruction efficiency between data and simulation [37].
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.

JHEP06(2021)044 6 Systematic uncertainties
The efficiency ratio in eq. (5.1) 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 [37]. Alternative parametrisations [36] 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 as described in refs. [38,39] with systematic uncertainties included with log-normal uncertainties. Upper limits on the branching fractions are determined using the observed distribution of CL s as a function of the assumed branching -11 -
A further 7.6 % uncertainty from the branching fraction of D + (s) → (φ → µ − µ + ) π + and a 1.0 % uncertainty from the finite size of the simulated 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.   fraction. The upper limits at 90 % and 95 % confidence level (CL) are given in table 3 and are shown in figure 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 and the results are computed under the assumption that the decay particles are uniformly distributed across the phase space. 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.

LHCb
Observed Expected ±1σ, ±2σ BaBar CLEO LHCb Figure 6. Upper limits at 90 % confidence level on the D + (s) signal channels. The median (orange), ±1σ and ±2σ expected limits are shown as box plots and the observed limit is given by a black cross. The green line shows the previous world's best limit for each channel where the solid, dotted and dashed lines correspond to BaBar, CLEO and LHCb [17][18][19]. The expected limits are determined using pseudo-experiments to approximate the CL s distribution under the backgroundonly hypothosis.