Search for the Rare Decays B 0 s → e + e − and B 0 → e + e −

A search for the decays B 0 s → e þ e − and B 0 → e þ e − is performed using data collected with the LHCb experiment in proton-proton collisions at center-of-mass energies of 7, 8, and 13 TeV, corresponding to integrated luminosities of 1, 2, and 2 fb − 1 , respectively. No signal is observed. Assuming no contribution from B 0 → e þ e − decays, an upper limit on the branching fraction B ð B 0 s → e þ e − Þ < 9 . 4 ð 11 . 2 Þ × 10 − 9 is obtained at 90(95)% confidence level. If no B 0 s → e þ e − contribution is assumed, a limit of B ð B 0 → e þ e − Þ < 2 . 5 ð 3 . 0 Þ × 10 − 9 is determined at 90(95)% confidence level. These upper limits are more than one order of magnitude lower than the previous values.

In this Letter, a search for B 0 s → e þ e − and B 0 → e þ e − decays is presented using data collected with the LHCb experiment in proton-proton collisions at center-of-mass energies of 7 TeV in 2011, 8 TeV in 2012and 13 TeV in 2015 to integrated luminosities of 1, 2, and 2 fb −1 , respectively. The signal yields are determined from a fit to the data and normalized to those of the B þ → J=ψK þ decay, where the J=ψ meson decays to e þ e − , which has a precisely measured branching fraction [12] and a similar dielectron signature in the detector.
The LHCb detector [13,14] 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. 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.
The online event selection is performed by a trigger [15], 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 high-energy deposit in the calorimeters associated with a signal electron candidate, or a muon candidate with high transverse momentum p T , or a photon, electron, or hadron candidate with high transverse energy from the decays of other particles from the pp collision. The software trigger requires a two-track secondary vertex with a significant displacement from any primary pp interaction vertex (PV). At least one charged particle must have high p T and be inconsistent with originating from a PV. A multivariate algorithm [16,17] is used in the trigger for the identification of secondary vertices consistent with the decay of a b hadron. Simulated samples are used to optimize the candidate selection, estimate selection efficiencies and describe the expected invariant-mass shapes of the signal candidates and background decays. In the simulation, pp collisions are generated using PYTHIA [18] with a specific LHCb configuration [19]. Decays of unstable particles are described by EVTGEN [20], in which final-state radiation is generated using PHOTOS [21]. The interaction of the generated particles with the detector, and its response, are implemented using the GEANT4 toolkit [22] as described in Ref. [23]. The simulation is corrected for data-simulation differences in B-meson production kinematics, detector occupancy, and isolation criteria [24] using B þ → J=ψK þ and B 0 s → J=ψϕ decays, with J=ψ → e þ e − and ϕ → K þ K − . Particle identification variables are calibrated using data from B þ → J=ψK þ and D 0 → K − π þ decays [25]. The calibration data are binned in momentum and pseudorapidity of the particle as well as detector occupancy to account for possible differences in kinematics between the investigated decay and the calibration data.
The B 0 ðsÞ → e þ e − candidates are selected in events passing the trigger requirements by combining two tracks that are inconsistent with originating from any PV in the event and which form a good-quality secondary vertex. The tracks are also required to have a momentum larger than 3 GeV=c and p T greater than 500 MeV=c , and must be identified as electrons using information from the Cherenkov detectors and calorimeters. The dielectron candidate's momentum must be aligned with the vector pointing from a PV (the associated PV) to the two-track vertex and have a considerable transverse component. The candidate must also have an invariant mass in the range ½4166; 6566 MeV=c 2 .
The measured electron momenta are corrected for losses due to bremsstrahlung radiation by adding the momentum of photons consistent with being emitted upstream of the magnet [26]. Candidates in data and simulation are separated into three categories with either zero, one, or both electrons having a bremsstrahlung correction applied. To avoid experimenters' bias, the narrowest dielectron invariant-mass region containing 90% of simulated B 0 s → e þ e − decays, corresponding to a range of ½4689; 5588 MeV=c 2 , was removed from the data set until the analysis procedure was finalized.
Candidates for the normalization mode, B þ → J=ψK þ , are constructed similarly, but require an additional track consistent with being a kaon and originating from the same vertex as the dielectron candidate. The dielectron candidate must have an invariant mass in the range ½2450; 3176 MeV=c 2 , consistent with arising from a J=ψ meson decay. In addition, the reconstructed B þ -candidate mass, when the dielectron candidate is constrained to the known J=ψ mass [12], must be above 5175 MeV=c 2 , suppressing partially reconstructed decays.
A boosted decision tree (BDT) algorithm [27][28][29] is used to separate B 0 ðsÞ → e þ e − signal from random combinations of two electrons (combinatorial background). The BDT is trained separately for data taking periods 2011-2012 (Run 1) and 2015-2016 (Run 2) on simulated B 0 s →e þ e − decays as signal proxy and dielectron candidates from data with a mass above 5588 MeV=c 2 as background proxy. The split between the data taking periods is done to account for changes in the center-of-mass energies and trigger strategies, which significantly impact the data distributions and improve the BDT and the particle identification algorithms in Run 2. It is checked that the data behave consistently across the data-taking periods. The BDT input variables comprise of the following: kinematic information on the electron tracks and B candidate, information on the displacement of the electrons and B candidate from the associated PV, and isolation variables that quantify the compatibility of other tracks in the event with originating from the same decay as the B candidate [24,30]. Candidates with a BDT response compatible with that of the background are discarded, with the threshold chosen by maximizing the figure of merit ϵ signal =ð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N background p þ3=2Þ [31], where ϵ signal is the signal efficiency and the expected background yield in the signal region is N background .
The final selected data set is separated by data-taking period and by category of bremsstrahlung correction. The branching fraction BðB 0 ðsÞ → e þ e − Þ is measured relative to that of the normalization channel via εðB 0 ðsÞ →e þ e − Þ and εðB þ →J=ψK þ Þ denote the efficiencies of the signal and normalization modes, and NðB 0 ðsÞ →e þ e − Þ and NðB þ → J=ψK þ Þ their yields. The normalization mode branching fraction (including that for the decay J=ψ → e þ e − )i sBðB þ → J=ψK þ Þ¼ð6.03 AE 0.17Þ × 10 −5 , taken from Ref. [12]. The b-hadron fragmentation fraction ratio f d =f u is assumed to be unity, while f s =f u ¼ 0.259 AE 0.015 [32] is used for the Run 1 data and is scaled by 1.068 AE 0.016 for the Run 2 data, according to Ref. [33], to account for center-of-mass energy differences. A measurement of f s =f u from Run 2 yields a consistent, but less precise, result [34].
The yield of the normalization mode is determined using an unbinned maximum-likelihood fit to the K þ e þ e − invariant mass separately for each year of data PHYSICAL REVIEW LETTERS 124, 211802 (2020) taking and bremsstrahlung category. The fit model comprises a Gaussian function with power-law tails [35] for the signal component, where the tail parameters are fixed from simulation, and an exponential function to describe combinatorial background. Summed over the bremsstrahlung categories, the yield of the normalization mode is 20480 AE 140 in the Run 1 data and 33080 AE 180 in the Run 2 data.
The selection efficiencies εðB 0 ðsÞ → e þ e − Þ and εðB þ → J=ψK þ Þ are determined separately for each year of data taking and bremsstrahlung category using simulated decays that are weighted to better represent the data. Calibration data are used to evaluate particle-identification efficiencies [25]. Trigger efficiencies are also estimated from data, using the technique described in Ref. [36].F o r simulated B 0 s → e þ e − decays, the mean B 0 s lifetime [37] is assumed. The selection efficiency is assumed to be the same for both B 0 → e þ e − and B 0 s → e þ e − decays, which is consistent with results from simulation. The normalization factors, α, are combined across the data-taking periods and, given in Table I, split by bremsstrahlung category (for the selection efficiency ratio between normalization and signal mode, see the Supplemental Material [38]).
In addition to the combinatorial background, backgrounds due to misidentification and partial reconstruction are present in the data. These backgrounds differ significantly between the categories of bremsstrahlung correction. Their invariant-mass shapes and relative contributions are evaluated using simulation. In the lower mass region, partially reconstructed backgrounds of the types B → Xe þ e − and B þ → D 0 ð→ Y þ e −ν e Þe þ ν e dominate, where X and Y represent hadronic systems. The main source of background in the B -mass region, however, stems from misidentified particles in the decays B 0 → π − e þ ν e and B → h þ h 0− , where h and h 0 are hadrons. The latter has a peaking structure in the B -mass region. Backgrounds involving misidentified particles contribute mostly to categories in which at most one of the electrons has a bremsstrahlung correction applied. The contribution from combinatorial background is evaluated from samesign lepton pairs in data and found to be small. The yields of the backgrounds are Gaussian constrained to their expected values, estimated from simulation using their known branching fractions [12].
The shape of the invariant mass of the B 0 s → e þ e − and B 0 → e þ e − components is modeled using a Gaussian function with power-law tails, where the parameters are obtained from simulation and differ between each bremsstrahlung category and year of data taking. The peak values and the widths of the functions are corrected for datasimulation differences by a factor determined from the normalization mode. The parameters of the B 0 s → e þ e − and B 0 → e þ e − line shapes are fixed to the same values with the exception of the peak value, which is shifted according to the known B 0 s -B 0 mass difference [12]. Due to the limited mass resolution, arising from imperfect bremsstrahlung recovery, the line shapes from B 0 s → e þ e − and B 0 → e þ e − are highly overlapping. Therefore the branching fraction of B 0 s → e þ e − is obtained by performing a simultaneous fit to the dielectron invariant-mass distribution of all six data sets while neglecting the contribution from B 0 → e þ e − , and vice versa. In these fits, the only shared parameters between categories are the branching fractions BðB 0 ðsÞ → e þ e − Þ and BðB þ → J=ψK þ Þ, and the ratio of the fragmentation fractions f s =f u .
Systematic uncertainties are estimated separately for each data set. Dominant sources of systematic uncertainties in the normalization arise from the uncertainty on the fragmentation fraction ratio, the technique used to evaluate the trigger efficiencies, and the determination of particleidentification efficiencies; the systematic uncertainties from these sources extend to 5.8%, 5.3%, and 5.3% on the branching fractions, respectively. The uncertainty on BðB þ → J=ψK þ Þ of 2.8% [12] is taken into account. A difference of up to 4.1% is found between the efficiency of the BDT selection on simulated B þ → J=ψK þ decays and B þ → J=ψK þ decays in data, which is assigned as a systematic uncertainty. The fraction of candidates in each bremsstrahlung-correction category of the signal modes is taken from simulation. The difference between simulation and data is investigated using B þ → J=ψK þ decays and its effect on the normalization, up to 4.0%, is taken as a systematic uncertainty. Systematic uncertainties on the invariant-mass resolution corrections are determined by repeating the correction procedure with pseudoexperiments obtained with the bootstrapping method [39], yielding up to 1.1%. A difference between the total selection efficiencies in the B 0 s → e þ e − and B 0 → e þ e − channels of up to 2.5% is assigned as a systematic uncertainty on the B 0 → e þ e − normalization factor. Due to the presence of an additional kaon in the final state of the normalization mode, the trackreconstruction efficiency is different between the signal and normalization modes. An uncertainty of 1.1% is assigned to the branching fraction as a systematic uncertainty on the kaon reconstruction efficiency arising from the limited knowledge of the interactions in the detector material [40]. Finally, an uncertainty of 1.0% is assigned to account for small differences in detector occupancy between the signal  (2020) 211802-3 and normalization mode arising from the trigger selection.
The dominant sources of systematic uncertainties on the background composition are due to the imprecise knowledge of the branching fractions of the background components. The largest uncertainty of this type on the expected background yield in the B-mass region is 14%, determined from refitting the mass sidebands while varying the background components according to their uncertainties. Taking all correlations into account, overall single event sensitivities of ½4.71 AE 0.12ðstatÞAE0.33ðsystÞ × 10 −10 for B 0 s →e þ e − and ½1.271AE0.034ðstatÞAE0.063ðsystÞ×10 −10 for B 0 → e þ e − are obtained.
The dielectron invariant-mass spectrum, summed over bremsstrahlung categories, is shown in Fig. 1, with  Upper limits on the branching fractions are set using the CL s method [41], as implemented in the GAMMACOMBO framework [42,43] with a one-sided profile likelihood ratio [44] as test statistic. The likelihoods are computed from fits to the invariant-mass distributions. In the fits, the normalization factor, normalization mode branching fraction, fragmentation fraction ratio, and background yields are Gaussian constrained to their expected values within statistical and systematic uncertainties. Pseudoexperiments, in which the nuisance parameters are set to their fitted values from data, are used for the evaluation of the test statistic.
In conclusion, a search for the rare decays B 0 ðsÞ → e þ e − is performed using data from proton-proton collisions  recorded with the LHCb experiment, corresponding to a total integrated luminosity of 5 fb −1 . No excess of events is observed over the background. The resulting limits on the branching fractions are BðB 0 s → e þ e − Þ < 9.4ð11.2Þ × 10 −9 and BðB 0 → e þ e − Þ < 2.5ð3.0Þ × 10 −9 at 90(95)% confidence level, when neglecting the contribution from the other decay. The mean B 0 s lifetime is assumed for B 0 s → e þ e − decays. Assuming SM-like CP-odd (CP-even) B 0 s → e þ e − decays, an increase (decrease) of 2.4% with respect to the quoted limit is found. The results improve the limits on these branching fractions [11] by more than one order of magnitude and constrain contributions beyond the SM, for example from scalar and pseudoscalar currents [10].
We express our gratitude to our colleagues in the CERN accelerator departments for the excellent performance of the LHC. We thank the technical and administrative staff at the LHCb institutes. We acknowledge support from CERN and from the national agencies: CAPES, CNPq,