Search for the decay $D_s^+\rightarrow \gamma e^+\nu_e$

A search for the rare radiative leptonic decay $D_s^+\to\gamma e^+\nu_e$ is performed for the first time using electron-positron collision data corresponding to an integrated luminosity of 3.19 fb$^{-1}$, collected with the BESIII detector at a center-of-mass energy of 4.178 GeV. No evidence for the $D_s^+\to\gamma e^+\nu_e$ decay is seen and an upper limit of $\mathcal B(D_s^+\to\gamma e^+\nu_e)<1.3\times 10^{-4}$ is set on the partial branching fraction at a 90\% confidence level for radiative photon energies $E_{\gamma}^*>0.01$~GeV.

A search for the rare radiative leptonic decay D + s → γe + νe is performed for the first time using electron-positron collision data corresponding to an integrated luminosity of 3.19 fb −1 , collected with the BESIII detector at a center-of-mass energy of 4.178 GeV. No evidence for the D + s → γe + νe decay is seen and an upper limit of B(D + s → γe + νe) < 1.2 × 10 −4 is set on the partial branching fraction at 90% confidence level for radiative photon energies E * γ > 0.01 GeV.

I. INTRODUCTION
In the Standard Model, the purely leptonic decays of heavy pseudoscalar mesons, P → e + ν e , are helicitysuppressed by a factor m 2 e . The helicity suppression in these processes can be overcome by the emission of a radiative photon as shown in Fig. 1. As a result, the decay rate of the purely leptonic radiative decay P → γe + ν e may be 10 3 − 10 5 times [1] larger than that of P → e + ν e . For example, the BFs of D + (s) → γe + ν e are theoretically predicted to range from 10 −5 to 10 −3 [2][3][4][5][6][7][8]. An experimental search for these decays can shed light on the dynamics of the underlying processes and can provide input of decay rates to theoretical calculations.
Previously, the BESIII experiment has searched for the radiative leptonic decay D + → γe + ν e using a data sample collected at a center-of-mass energy √ s = 3.773 GeV. No significant signal is observed and an upper limit on the partial decay BF for radiative photon energies E * γ > 0.01 GeV is set to B < 3.0 × 10 −5 at the 90% confidence level (C.L.) [9], approaching the range of theoretical predictions, (1.9 − 2.8) × 10 −5 [5,6]. The decay D + → γe + ν e is Cabibbo-suppressed, while the decay D + s → γe + ν e is Cabibbo-favored. The full BF of D + s → γe + ν e is predicted to be of the order 10 −5 − 10 −4 in the light front quark model [2] and in the nonrelativistic constituent quark model [4]. The theoretical study in Ref. [5] indicates that the long-distance contribution described by the vector meson dominance model, as shown in Fig. 2, may further enhance this decay BF up to order 10 −4 . Moreover the BF is predicted to be of order 10 −3 within the perturbative quantum chromodynamics method combining heavy quark effective theory [3]. With a BF sensitivity of 10 −4 − 10 −5 , this decay may be detectable at BESIII.
In this paper, we report on the first search for the radiative leptonic decay D + s → γe + ν e , using a data sample corresponding to an integrated luminosity of 3. 19  analysis procedure of the nominal analysis has been developed as a blind analysis, based on an inclusive Monte Carlo (MC) simulated data sample with equivalent luminosity the same as data. The inclusion of the charge conjugate process is implied throughout the paper unless explicitly specified otherwise.

II. BESIII DETECTOR AND DATA SET
The BESIII detector is a magnetic spectrometer [10] located at the Beijing Electron Positron Collider (BEPCII) [11]. The cylindrical core of the BESIII detector consists of a helium-based multi-layer drift chamber (MDC), a plastic scintillator time-of-flight system (TOF), and a CsI(Tl) electromagnetic calorimeter (EMC), which are all enclosed in a superconducting solenoidal magnet providing a 1.0 T magnetic field. The solenoid is supported by an octagonal flux-return yoke with resistive plate counter muon identifier modules interleaved with steel. The acceptance of charged particles and photons is 93% over 4π solid angle. The chargedparticle momentum resolution at 1 GeV/c is 0.5%, and the specific energy loss (dE/dx) resolution is 6% for the electrons from Bhabha scattering. The EMC measures photon energies with a resolution of 2.5% (5%) at 1 GeV in the barrel (endcap) region. The time resolution of the TOF barrel part is 68 ps. The end cap TOF system was upgraded in 2015 with multi-gap resistive plate chamber technology, providing a time resolution of 60 ps [12,13].
MC simulated events are generated with the geant4based [14] software package boost [15] that describes the detector geometry and material, implements the detector response, simulates digitization, and incorporates time-dependent beam backgrounds. An inclusive simulation sample, which includes open charm processes, the initial state radiation (ISR) production of ψ(3770), ψ(3686) and J/ψ, qq (q = u, d, s) continuum processes, along with Bhabha scattering and µ + µ − , τ + τ − and γγ processes, is produced at √ s = 4.178 GeV. The open charm processes are simulated using conexc [16]. The effects of ISR and final state radiation (FSR) [17] are taken into account. Decays of unstable particles are simulated by evtgen [18] using BFs from the Particle Data Group [19], and the remaining unknown decay modes of ψ are generated using the modified lund model [20]. The signal candidates are simulated using the method employed in Ref. [9], where the two parameters, the decay constant [19] and the quark mixing matrix element [19] are adjusted according to the decay channel. The minimum energy of the radiative photon is set at 0.01 GeV to avoid the background of soft photons caused by infrared divergence.

III. DATA ANALYSIS
At √ s = 4.178 GeV, the D s mesons are mostly produced in the process e + e − → D + s D * − s . This allows to perform the analysis using a modified double-tag (DT) technique [21]. First, the D − s decay is fully reconstructed, leading to the single-tag (ST) events. In these ST events, the signal decay D + s → γe + ν e , called the DT, is selected and investigated in the presence of one additional isolated photon or π 0 meson originating from the D * − s decay. The BF of D + s → γe + ν e is determined by where N tot ST and N signal are the ST and DT yields in data, respectively. ǫ γ soft (π 0 soft )SL is the reconstruction ef- soft ) denotes the soft γ or π 0 from the D * − s , γe + ν e decays come from either the bachelor D + s or D * + s , ǫ i ST and ǫ i DT are the efficiencies of selecting the ST and DT candidates, and i denotes the i-th tag mode as described below.
The ST candidates are reconstructed through the decay modes D − where the subscripts of η (′) represent the decay modes used to reconstruct η (′) . All charged tracks must have a polar angle (θ) within | cos θ| < 0.93. The reconstructed tracks are required to point back to the interaction point (IP) region with |V r | < 1 cm and |V z | < 10 cm, where |V r | and |V z | are the distances of closest approach to the IP in the transverse plane and along the positron beam direction, respectively. Charged kaons and pions are identified by using the combined information from dE/dx and TOF. The charged tracks are assigned as pion (kaon) candidates if L π(K) > L K(π) , where L π(K) is the confidence level for the pion (kaon) hypothesis. Below 1.2 GeV/c, the particle identification (PID) efficiencies of charged kaons (pions) range from 89% (85%) to 99%, while the rates of misidentifying kaons (pions) as pions (kaons) range from 1% to 12% (15%).
The K 0 S candidates are formed from pairs of oppositely charged tracks satisfying |V z | < 20 cm. The two charged tracks are taken as π + π − without identification requirements and are constrained to have a common vertex. The invariant mass M of the π + π − pair is required to be within (0.487, 0.511) GeV/c 2 . The decay length of the K 0 S candidate is required to be larger than twice the vertex resolution away from the IP.
Photon candidates are reconstructed from clusters of energy deposited in the EMC, with the energy measured in nearby TOF counters included to improve reconstruction efficiency and energy resolution The energies of photon candidates must be larger than 0.025 GeV for the barrel or 0.05 GeV for the endcap region. These requirements are safe for the minimum energy requirement E * γ > 0.01 GeV on the radiative photon. The cluster timing [22] is required to be between 0 and 700 ns to suppress electronic noise and energy depositions unrelated to the event of interest.
To remove the soft pions coming from D * decay, the momentum of the pion coming directly from the ST D − s decay must be larger than 0.1 GeV/c. For the π + π − π − and K − π + π − final states, the contributions [19]. The ST D − s mesons are identified by the modified mass and the D − s recoil mass where p D − s is the three-momentum of the ST candidate in the rest frame of the e + e − system, M D − s is the nominal D − s meson mass [19] and E beam is the beam energy. The non-D + s D * − s events are suppressed by requiring M mod ∈ (2.010, 2.073) GeV/c 2 . In each event, only the candidate with the M rec closest to the D * + s nominal mass [19] is chosen. The invariant mass (M tag ) spectra of the accepted ST candidates for the 14 tag modes are shown in Fig. 3. The ST yield is determined via unbinned maximum-likelihood fits to each spectrum. Signals and the D − → K 0 S π − peaking background with a tiny fraction (dashed-black line in Fig. 3) in the D − s → K 0 S K − mode are described by MC-simulated shapes using the kernel density estimation method [23]. To take into account the resolution difference between data and simulation, the MC-simulated shapes are convolved with a Gaussian function for each tag mode, where the parameters of the Gaussian function are left free in the fit. The non-peaking background is modeled by a secondor third-order Chebychev polynomial function, and the reliability of the fitted non-peaking background has been verified using the inclusive MC sample. Candidates in the signal regions, denoted by the boundaries in each sub-figure of Fig. 3, are kept for further analysis. The M tag signal regions, the ST yields in data and the ST efficiencies are summarized in Table I. The total ST yield is N tot ST = 395412 ± 1931, where the uncertainty is statistical.
The D + s → γe + ν e candidates are selected from the remaining charged tracks and showers in the side recoiling against the ST D − s meson and the isolated photon or π 0 meson with the same criteria as used in the ST candidate selection. It is required that there is only one good charged track, with charge opposite to the ST D − s meson. The positron is identified using the confidence level computed by combining PID information from dE/dx, TOF and EMC. Under the assumption that the charged track in the signal decay is a positron, a pion or a kaon, three confidence levels are calculated: L ′ e , L ′ π and L ′ K . The charged track is identified as a positron if L ′ e > 0.001 and L ′ e /(L ′ e + L ′ π + L ′ K ) > 0.8. To reduce the rate of misidentifying a pion as a positron, the ratio E e /p e is required to be greater than 0.8, where E e and p e are the deposited energy of the positron in the EMC and the momentum measured by the MDC, respectively. Below 1.2 GeV/c, the PID efficiencies of e ± are greater than 98%, while the averaged rate of misidentifying K ± or π ± as e ± is about 0.3%. To improve the degraded momentum resolution of the electron due to FSR and bremsstrahlung effects, the energies of neighboring photons are added back to the positron candidates. Specifically, the photons with energy greater than 0.03 GeV and within a cone of 5 • around the positron direction (but excluding the radiative photon candidate) are included.
To select the radiative semileptonic decay candidate from the process e + e − → D + candidates are constrained to their individual nominal masses. In addition, the neutrino is treated as a missing particle in the DT event. The hypothesis with the smallest χ 2 kine is chosen. The χ 2 kine distribution of the accepted candidates is shown in Fig. 4.
To suppress the background from D + s hadronic decays due to fake photons and charged tracks, the maximum energy of the showers not used in the DT event selection (E max γ extra ) is required to be less than 0.2 GeV, and events with additional charged tracks (N extra char ) are removed. To suppress backgrounds from D + s → τ + ν τ and D + s → ηe + ν e , χ 2 kine is required to be less than 70. The backgrounds from D + s → ηe + ν e are further suppressed by rejecting the events if the invariant mass of any γγ combination that has not been used in ST selection satisfies M γγ ∈(0.51, 0.56) GeV/c 2 . These requirements keep 80% of the signal events, but remove more than 70% of the background events.
Finally, the signal candidates are searched for in the data distribution of the kinematic variable where p miss ≡ −( p γ + p e + p tag + p γ soft (π 0 soft ) ) (6) in the e + e − rest frame. The distribution of U miss of the surviving DT candidates is shown in Fig. 5. The signal Mode  candidates of D + s → γe + ν e should peak around zero in the U miss distribution, as shown by the signal MC sample (black dashed line). Figure 6 shows the E γ distribution in the U miss signal region (−0.06, 0.06) GeV, where the data points overlap with the simulated distributions of the backgrounds coming from the D + s → ηe + ν e and D + s → τ + ν τ decays. No excess of signal candidates is observed in the region where they are expected to peak.

IV. RESULT
To measure the signal yield of the D + s → γe + ν e decay, an extended unbinned maximum likelihood fit is performed to the U miss distribution. The result of the fit is shown as the solid line in Fig. 5. The signal and background shapes of the decays D + s → ηe + ν e with η → γγ and D + s → τ + ν τ with τ + → e + ν eντ are determined from the corresponding MC sample for each decay mode. For the other background components, the shape is determined from the inclusive simulation sample. The DT efficiencies of the individual ST modes are listed in Table I. We obtain a signal yield of N signal = −14 ± 15. Since no significant signal is observed, an upper limit on the BF of the D + s → γe + ν e decay at the 90% C.L. is set by solving the equation [19] B UL 0 L(B)dB = 90%.
A series of fits on the U miss distribution is carried out, fixing the BF at different values. The resulting likelihood distribution L is shown in Fig. 7. The upper limit on the BF at the 90% C.L. is found to be 7.1 × 10 −5 .
The sources of systematic uncertainties that affect the upper limit calculation are discussed below. With the DT method, the systematic uncertainties related to the selection of the ST candidates are found to be negligible. To estimate the uncertainty in the ST yield and to avoid statistical fluctuations, a total of 1000 fits to generated samples have been performed by using alternative signal (double Gaussian function) and background (Chebyshev polynomial) shapes. The systematic uncertainties of 0.3% and 0.2% are obtained by taking the mean value of the distribution of the relative normalized difference between the pseudo-experiments and baseline fit results. The total systematic uncertainty in the ST tag yield is taken as the squared sum, and it is found to be 0.4%. To estimate the systematic uncertainty due to the D + s → γe + ν e form factors, an alternative signal MC sample based on the single-pole model [6] has been pro-duced, the difference between the DT efficiency obtained with this model and the one with our nominal model at 0.025 GeV is 2.6%, and the difference of fractions of the generated events in (0.01, 0.025) GeV between two models is 8%. Due to full correlation of the two systematic errors, they are added linearly to obtain the systematic uncertainty in form factor model, 11%. The systematic uncertainties attributed to the positron tracking and PID efficiencies are studied with a control sample of radiative Bhabha scattering events. The control sample and the D + s → γe + ν e simulation sample have different distributions in the momentum and angle of the positron. To account for these differences, a correction resulting from a two-dimensional re-weighting in momentum and angle is applied to the positron tracking efficiency and to the positron PID efficiency. The total systematic error caused by uncertainties in positron tracking and PID is estimated to be 0.4%. The systematic uncertainty in the photon selection is evaluated using a control sample of J/ψ → π + π − π 0 decays [24]. It is determined to be 1.0%. Systematic uncertainties of 1.1% and 0.9% due to the E max γ extra and N extra char selection criteria are estimated by analyzing the DT hadronic D * + s D − s events. Systematic uncertainties of 2.2% and 0.3% due to the χ 2 kine requirement and the FSR effect are computed by repeating the fit without the χ 2 kine requirement and the correction for the FSR effect, respectively, and taking the differences with respect to the baseline fit. To estimate the uncertainty of U miss fitting related to the background shape, the fraction of each of the main background components is varied within one standard deviation of the corresponding BF [19]. The largest deviation with respect to the baseline result is 10%. To avoid statistical fluctuations, a study based on pseudo−experiments is performed. A total of 1000 fits to generated samples is performed by varying the background shape. A systematic uncertainty of 10% is obtained by taking the mean value of the distribution of the relative normalised difference between the pseudo−experiments and the baseline fit results. Differences between the tag of the ST modes in data and simulation are expected to impact the final result due to the different multiplicities. The associated systematic uncertainty is assigned as 0.5% by studying the tracking/PID efficiencies and the photon selection in different multiplicities resulting in a difference between data and MC sample. Table II summarizes all the systematic uncertainties. The impact of the systematic uncertainty on the upper limit of the BF is taken into account by convolving the distribution of the sensitivity (S) where LH(t) = Cexp −(t−t) 2 2σ 2 t , C is a normalization constant,t and σ t can be obtained when the likelihood distribution is fitted by LH(t). The valueŜ is the nominal efficiency and δ S is the systematic uncertainty on the BF [25]. Finally, the upper limit on the BF of the D + s → γe + ν e decay is set to be 1.2 × 10 −4 at the 90% C.L..

V. SUMMARY
The first search for the radiative leptonic decay D + s → γe + ν e is performed using e + e − collision data corresponding to an integrated luminosity of 3.19 fb −1 collected at √ s = 4.178 GeV, by employing a DT technique. No significant signal for the signal decay D + s → γe + ν e is observed. With a 0.01 GeV cutoff on the radiative photon energy, the upper limit on the BF of the D + s → γe + ν e decay mode is set to be B(D + s → γe + ν e ) < 1.2 × 10 −4 at