Search for the Higgs boson decay to a pair of electrons in proton-proton collisions at

A search is presented for the Higgs boson decay to a pair of electrons ( e + e − ) in proton-proton collisions at √ s = 13 TeV. The data set was collected with the CMS experiment at the LHC between 2016 and 2018, corresponding to an integrated luminosity of 138fb − 1 . The analysis uses event categories targeting Higgs boson production via gluon fusion and vector boson fusion. The observed upper limit on the Higgs boson branching fraction to an electron pair is 3 . 0 × 10 − 4 (3 . 0 × 10 − 4 expected) at the 95% conﬁdence level, which is the most stringent limit on this branching fraction to date. © 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons .org /licenses /by /4 .0/). Funded by SCOAP 3 .


Introduction
Since the discovery of the Higgs boson by the ATLAS and CMS Collaborations in 2012 [1][2][3], various measurements of its interactions with standard model (SM) particles have been performed.In particular, the interactions of the Higgs boson with the thirdgeneration fermions have been observed, with coupling strengths found to be in agreement with the expectations of the SM [4][5][6][7][8][9].More recently, the CMS Collaboration has also presented first evidence for the decay of the Higgs boson to the second-generation fermions, in a search for H → μ + μ − decays [10].Such measurements of Yukawa couplings to light fermions are experimentally challenging at the LHC, given that the SM predicts the coupling strengths to be proportional to the fermion mass [11][12][13].Consequently, couplings of the Higgs boson to the first-generation fermions have yet to be confirmed experimentally.
For the Higgs boson decay to an electron pair, H → e + e − , the SM predicted branching fraction of B(H → e + e − ) ≈ 5 × 10 −9 [14] is extremely small and inaccessible with the current LHC data set.Higher-order corrections included in the branching fraction prediction introduce diagrams which do not contain the Higgs boson Yukawa coupling to the electron.These contributions are expected to dominate in the SM, but are still much smaller than present sensitivity.Therefore, a search for this decay channel currently provides the only direct probe of the Higgs boson Yukawa coupling to the electron, which is enhanced in several scenarios beyond the E-mail address: cms -publication -committee -chair @cern .ch.SM.The simplest of these are models postulating two Higgs doublets [15]; other extensions include the addition of higher order operators to the SM Lagrangian, including dimension 10 operators that could modify the coupling by a factor of ≈10 [16,17].These enhancements are, however, still smaller than present sensitivity could detect.This Letter describes a search for H → e + e − decays in protonproton (pp) collisions at √ s = 13 TeV, using data recorded with the CMS detector in 2016-2018, corresponding to an integrated luminosity of 138 fb −1 .The most recent search from the CMS Collaboration for this decay was performed using pp collision data collected at √ s = 8 TeV with an integrated luminosity of 19.7 fb −1 [18].A 95% confidence level (CL) upper limit on B(H → e + e − ) of 1.9 × 10 −3 was determined.The most sensitive search for H → e + e − decays from the ATLAS Collaboration was performed using pp data with an integrated luminosity of 139 fb −1 , collected at √ s = 13 TeV.A 95% CL upper limit on B(H → e + e − ) was set at 3.6 × 10 −4 [19].
The analysis presented in this Letter primarily targets Higgs boson production by gluon fusion (ggH) and vector boson fusion (VBF).Rarer Higgs boson production processes, including production in association with top quarks (ttH) and a vector boson (VH), are also considered, although no dedicated selection is applied to target them.The final states of interest comprise a pair of two prompt, oppositely charged electrons, potentially produced in association with hadronic jets.The sensitivity to each production mode is increased by dividing events into analysis categories, giving an improved signal to background ratio (S/B).Each analysis category is defined by a selection on a dedicated multivariate (MVA) dis-criminant.The MVA-based classifiers are trained to separate Higgs boson signal events from the dominant Drell-Yan (DY) background and smaller background contributions from tt decays.A simultaneous fit to the dielectron invariant mass distribution (m ee ) in each category is then used to extract an upper limit on B(H → e + e − ).

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections.The pseudorapidity (η) coverage provided by the barrel and endcap detectors is extended by the forward hadron calorimeter, which uses steel as an absorber and quartz fibres as the sensitive material.Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
Events of interest are selected using a two-tiered trigger system [20,21].The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed latency of about 4 μs.The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage.
The global event reconstruction (also called particle-flow (PF) event reconstruction [22]) aims to reconstruct individual particles (photons, charged and neutral hadrons, muons, and electrons) by optimally combining information from the various elements of the CMS detector.In the case of electrons, the energy is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the charged electron track.The momentum resolution for electrons with transverse momentum p T ≈ 45 GeV from Z → e + e − decays ranges from 1.6-5.0%[23,24].The candidate vertex with the largest value of summed physics-object p 2  T is taken to be the primary pp interaction vertex.The physics objects used for this determination are the jets, clustered using the infrared-and collinear-safe anti-k T algorithm [25,26] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector p T sum of those jets.
For charged hadrons, the energy is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for the response function of the calorimeters to hadronic showers.The energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.Hadronic jets are built from PF candidates using the anti-k T algorithm [25] with a distance parameter of 0.4.The momentum of a jet is determined as the vectorial sum of all particle momenta in the jet, and is found from simulation to be, on average, within 5-10% of the true momentum over the whole p T spectrum and detector acceptance.Charged hadrons originating from additional pp interactions are removed from the analysis.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [27].

Analysis strategy
The analysis strategy is similar to that used in the SM H → γγ analysis [28].In this case, however, events are selected using a dielectron trigger with asymmetric electron p T thresholds of 23 and 12 GeV, and loose isolation requirements.Offline, all events entering the analysis are required to pass a loose selection placed on the electron kinematics and an electron identification (ID) criterion.Simulated signal and background events are then used to train MVA-based classifiers designed to reduce the number of background events entering the selection.This is performed independently for the ggH and VBF channels, with classifiers trained on simulated samples corresponding to the three data-taking years, combined in accordance with their integrated luminosities.
In both ggH and VBF channels, the dominant background consists of DY events, with smaller contributions from tt processes with dilepton or lepton + jets final states.In the VBF channel, electroweak Z boson production in association with two jets (EW Z → e + e − ) is also considered.These events closely mimic the VBF signal process and are included as a background during MVA training.Although the cross section for VBF production is roughly an order of magnitude smaller than for ggH production, the VBF final state provides a distinct signature in the detector, in which the Higgs boson is produced in association with two forward jets with a large η difference and high dijet invariant mass.These features are exploited by the VBF classifier to enhance the expected S/B ratio, allowing the sensitivity of the VBF categories to become comparable to those targeting ggH.
The output scores of the classifiers are used to define analysis categories.The number and location of category boundaries is chosen to optimize the expected significance of all resulting categories combined.The categories are non-overlapping by construction, with events first being considered for the VBF categories, as these have both a higher S/B and a distinctive topology.If an event fails the VBF requirements, it is then considered for the ggH analysis categories.
Once the analysis categories are defined, a maximum likelihood fit is performed to the m ee distribution in each category simultaneously.The smoothly falling background spectrum, consisting mostly of DY events, is modelled directly from data using the discrete profiling method [29].The signal model is derived from simulated Higgs boson events independently for each production mode and analysis category.The fit is used to extract 95% CL upper limits on B(H → e + e − ).

Data and simulated events
This analysis uses pp collision data corresponding to an integrated luminosity of 138 fb −1 , of which 36.3, 41.5, and 59.8 fb −1 were collected during 2016, 2017, and 2018, respectively.
A Monte Carlo (MC) simulated sample is generated for each of the ggH, VBF, VH, and ttH signal processes using Mad-Graph5_amc@nlo v2.6.5 [30] at next-to-leading order (NLO) accuracy in perturbative quantum chromodynamics (QCD).For the ggH process, up to three additional partons are considered in the matrix element calculation; for VBF production, up to one additional parton is considered, excluding the quark-initiated jets produced in the leading order process.In addition, events from ggH production are weighted to match the predictions of the nnlops program [31], as a function of the Higgs boson transverse momentum and the number of jets in the event.Simulated DY and electroweak Z → e + e − background processes are also generated at NLO using MadGraph5_amc@nlo, while tt events are produced using the powheg v2.0 [32][33][34][35] NLO generator.All parton-level simulated samples are interfaced with pythia v8.230 [36] for parton showering and hadronization, using the CP5 underlying event tune [37].Parton distribution functions (PDFs) are taken from the NNPDF 3.1 [38] set for the simulation of all years.The ggH, VBF, VH, and ttH production cross sections used to normalize the signal samples are taken from the LHC Higgs Working Group recommendations [39], at the highest available order.Additional pp The CMS Collaboration Physics Letters B 846 (2023) 137783 interactions occurring in both the same and adjacent bunch crossings (pileup) are simulated as a set of minimum bias interactions that are mixed with the hard scattering event.In the analysis of each data-taking year, the simulated events are weighted based on the number of pileup events to match the distribution measured in data.Finally, the response of the CMS detector is simulated using the Geant4 package [40].

Event selection
A loose preselection on the events is applied to ensure they are consistent with a Higgs boson decaying to two electrons.The two electrons with the highest and next-highest p T are referred to as the leading and subleading electrons, respectively, and define the dielectron system.Each event must contain at least two oppositely charged electrons fulfilling the following requirements: (25) GeV for the leading (subleading) electron; • a tight working point on an electron ID, designed to have a 90% signal selection efficiency.The ID combines descriptions of the shape and isolation of energy deposits in the ECAL, the quality of associated tracks, and the compatibility of measurements from the tracker and ECAL, in an MVA-based discriminator [24]; • pseudorapidity within the ECAL acceptance (|η| < 2.5), and not in the barrel-endcap transition region (1.44 < |η| < 1.57); and • 110 < m ee < 150 GeV, chosen to limit contributions from DY events, and to ensure the m ee sideband regions have sufficient events to constrain the background expectation in the signal region.
All jets entering the analysis are required to have p T > 25 GeV and |η| < 4.7.Jets are required to pass a tight pileup identification criterion [41] that uses the topology of the jet shape, the number of charged and neutral jet constituents, and information on any associated tracks, to reject jets resulting from pileup.A tight requirement is also placed on an additional jet ID, which rejects spurious jets resulting from detector noise [41].Finally, a selection is applied to suppress the observed noise in the ECAL endcaps for low-p T jets in the 2017 data; this selection vetoes jets with p T < 50 GeV within the region 2.7 < |η| < 3.1.
Scale factors are applied to simulated samples to correct for differences between simulation and data.These include corrections for differences in the electron reconstruction and identification efficiencies, as well as in the trigger efficiency.

Event categorization
The categorization targeting VBF events are based on a boosted decision tree (BDT) trained to discriminate between signal and background events, referred to as the VBF BDT.Prior to training the BDT, an additional selection is applied on top of the loose preselection in order to target VBF events.This is referred to as the VBF preselection and comprises the following: • p T > 40 (25) GeV for the leading (subleading) jet; • dijet invariant mass > 350 GeV.
All simulated VBF events passing the VBF preselection are considered as signal when training the BDT.The largest background process in the VBF phase space consists of DY events, with a smaller but significant contribution from tt production.Although small in yield, the electroweak Z boson production process is also considered, since such events typically have signal-like characteristics and thus are assigned high VBF BDT scores.
Many observables characterising the VBF event topology are leveraged in the BDT in order to significantly increase the S/B in each analysis category, including information related to the kinematics of the individual electrons, the dielectron system, the three leading jets individually, the dijet system, and the dielectron plus dijet system.Inputs with good discriminating power include the invariant mass of the dijet object, the difference in η between the two leading jets, and a centrality variable [42] defined as , where η ee is the dielectron pseudorapidity, and j 1(2) labels the leading (subleading) jet.Other jet variables such as a quark-gluon identification score [43] are also included.The score is based on quantities such as the p T and multiplicity of PF candidates reconstructed within a jet in a likelihood discriminant to separate gluon and quark-initiated jets.Finally, it was verified that information related to the dielectron system was not sufficient for the Higgs boson mass (m H ) to be inferred by the classifier.The output score for the VBF BDT is shown for both simulation and data in Fig. 1 (left).
In this analysis, the background model is taken directly from data, whereas the signal modelling uses simulated samples.Therefore, it is necessary to ensure there is reasonable agreement between data and simulation for input observables to the VBF classifier as well as the output score.This is checked using Z → e + e − events from a control region, defined to be non-overlapping with the analysis category phase space.The VBF preselection is also applied, with the exception of the m ee requirements which are shifted to approximately centre on the Z boson mass (80 < m ee < 100 GeV).Drell-Yan events are chosen for this validation since the final state mimics that of the H → e + e − process and is relatively free from contaminating backgrounds.The distribution of the classifier output score is shown in Fig. 1 (right).Agreement between data and simulation is observed to be within the uncertainty given by the combined systematic and statistical variations, for the range of output scores in which the analysis categories are constructed.
Analysis categories targeting VBF events are defined using the output score of the BDT.The analysis defines only two VBF categories, since the sensitivity to VBF events is observed to saturate at numbers larger than two.The expected signal and background yields for both VBF categories are shown in Table 1.
The categorization for ggH events is also based on a BDT trained to discriminate between simulated signal and background events, referred to as the ggH BDT.All simulated signal and background events passing the basic preselection are considered when training.The largest background contribution in this channel consists of DY events, the kinematics of which are typically similar to the ggH signal.Hence, the classification task is nontrivial and many features are leveraged to improve the separation power of the BDT.These features include kinematic properties of the individual electrons, the dielectron system, and up to two jets, if present.Although these inputs are broadly similar between simulated ggH signal and background events, quantities such as the dielectron p T provide reasonable discriminating power, where the average dielectron p T for ggH events is higher than in background.Similarly to the VBF BDT, it is checked that information related to the dielectron system is not sufficient for m H to be inferred by the classifier.The output score for the ggH BDT is shown for both simulation and data in Fig. 2

(left).
The validation procedure for the ggH BDT is identical to that for the VBF BDT, described above, and is performed in the Z boson mass control region designed to be non-overlapping with the ggH analysis categories.The basic analysis preselection is also applied, with the exception of the m ee requirements.The distribution of the classifier output score in this region is shown in Fig. 2 (right).The agreement between data and simulation is typically within the uncertainty permitted by the combined systematic and statistical variations; any residual differences are smaller than the theoreti-Fig.1. Distribution of the output score of the VBF BDT in all simulated background and signal events, and data (left), passing VBF preselection.The ggH and VBF signals are scaled for better visibility.Category boundaries targeting VBF production are denoted with dashed lines.The shaded region defines events which are not selected to enter VBF analysis categories, but may populate those targeting ggH.A small excess in data with respect to the MC background prediction is observed in the most signal-like bin; however, these events are not found to peak around the signal region of the m ee distribution.The right plot shows the distribution of the output score of the VBF BDT in a control region around the Z boson mass.The simulated samples include both DY and EW Z → e + e − events.The combined impact of the statistical and systematic uncertainties in simulation is shown by the red shaded band, where the systematic component includes uncertainties in the jet energy scale and resolution corrections, alongside the electron energy scale corrections.Uncertainties in the efficiency of electron identification, reconstruction, and trigger selection are also included, as well as the uncertainty on the integrated luminosity, presented in Section 7. Good agreement is observed between the DY simulation (filled histogram) and data (black markers), within the phase space in which the analysis categories are constructed.

Table 1
The total expected number of signal events for m H = 125.38GeV in analysis categories targeting ggH and VBF events, for an integrated luminosity of 138 fb −1 .The fractional contribution from each production mode to each category is also shown.The σ eff , defined as the smallest interval containing 68.3% of the m ee distribution, is listed for each analysis category.The final column shows the expected ratio of signal to background, where S and B are the numbers of expected signal and background events in a ±1σ eff window centred on m H .The expected background event yields are extrapolated from data sidebands.cal uncertainty in the ggH production cross section, presented in Section 7. Four analysis categories targeting ggH events are defined using the output score of the ggH BDT.It is checked that going beyond four categories offers no improvement in sensitivity.The expected signal and background yields in each ggH category are shown in Table 1.

Systematic uncertainties
In this analysis, the systematic uncertainty associated with the background estimation from data is handled using the discrete profiling method, as described in Section 8. Systematic uncertainties that affect the signal model are implemented in one of the following two ways.Uncertainties that modify the shape of the m ee distribution are incorporated into the signal model as nuisance parameters.These are typically experimental uncertainties related to the energy measurement of the individual electrons.Conversely, if the shape of the m ee distribution is unaffected, the uncertainty is treated as a log-normal constrained variation in the event yield.

Theoretical uncertainties
The sources of theoretical uncertainty considered in this analysis are listed below.The effects of theoretical uncertainties are taken to be correlated across years.
• Renormalization and factorization scale uncertainty: the uncertainty arising from variations of the renormalization (μ R ) and factorization (μ F ) scales used when computing the expected SM cross section and event kinematics.These account for the missing higher-order terms in perturbative calculations.The recommendations provided by the Higgs Cross Section Working Group quoted in Ref. [39] are followed.The uncertainties in the signal acceptance due to the μ R and μ F scales are estimated using three sources: varying the μ R scale by a factor of 2 and 0.5, varying the μ F scale by a factor of 2 and 0.5, and varying both in the same direction simultaneously.The impacts of the signal acceptance uncertainties are evaluated keeping the overall normalization of each signal process constant, and are at largest (for the ggH production mode) 10%.Events with scores in the grey shaded region are discarded from the analysis.The right plot shows the distribution of the output score of the ggH BDT in a control region around the Z boson mass.Agreement is compared between the DY simulation (filled histogram) and data (black points).The combined impact of the statistical and systematic uncertainties in simulation is shown by the red shaded band, where the sources contributing to the systematic component are identical to those included in Fig. 1.Residual differences between data and simulation are smaller than the theory uncertainties in the ggH cross section and event kinematics, which are included in the final maximum likelihood fit.
• PDF uncertainties: these account for the uncertainty due to imperfect knowledge of the composition of the proton, that affects which partons are most likely to initiate high energy events.The normalization uncertainties are computed following the PDF4LHC prescription [14,44].The impact on the normalization ranges between 1.9 and 3%.
• Uncertainty in the strong force coupling constant: the uncertainty in the value of the strong force coupling constant α S is included in the treatment of the PDF uncertainties, following the PDF4LHC [14,44] prescription.The impact on the normalization is largest for ggH production, with a value of 2.6%.
• Underlying event and parton shower uncertainties: these uncertainties are obtained using dedicated simulated samples which vary the pythia tune from that used in the nominal simulation samples, and vary the renormalization scale for QCD emissions in both initial-state and final-state radiation by a factor of 2 and 0.5.These uncertainties are treated as migrations of events from a given production mode into and out of the VBF and ggH analysis categories.The largest effect comes from the parton shower uncertainty in VBF events, which can change the signal yield in the VBF analysis categories by up to 5%.
The total impact on the B(H → e + e − ) measurement from theoretical systematic uncertainties is +0.2 −0.1 × 10 −4 , with dominant contributions from uncertainties affecting the ggH production cross section.

Experimental uncertainties
The uncertainties that affect the shape of the signal m ee distribution are listed below.
• Electron energy scale: the uncertainty associated with the corrections applied to the electron energy scale in simulated events [24], derived using a sample of Z → e + e − tag-andprobe (T&P) events [45].Four nuisance parameters are defined for the possible combinations of low R 9 , high R 9 , electrons in the endcap, and electrons in the barrel region of the ECAL, where R 9 is defined as the energy sum of the 3 × 3 crystals centred on the most energetic crystal in the candidate electromagnetic cluster divided by the sum of the crystal energies assigned to that electromagnetic cluster.
• Electron energy scale nonlinearity: a further source of uncertainty to cover possible differences between the linearity of the electron energy scale between data and simulation, estimated on a sample of Z → e + e − T&P events.An uncertainty of 0.1% is assigned for electrons with p T < 80 GeV, and 0.2% for those above.
The uncertainties that only modify the event yield include: • Integrated luminosity: the integrated luminosities for the 2016, 2017, and 2018 data-taking years have individual uncertainties of 1.2, 2.3, and 2.5%, respectively [46][47][48].The total uncertainty for the 2016-2018 period is 1.6%.These are partially correlated across the different data sets to account for common sources of uncertainty in the luminosity measurement schemes.
• Electron identification and reconstruction: uncertainties on the scale factors derived to correct for differences in simulation and data for the electron identification and reconstruction efficiency.For both sources, the size of the uncertainty is approximately 1% in each category.
• Jet energy scale and resolution corrections: the energy scale of jets is measured using the p T balance of jets with Z bosons and photons in Z → e + e − , Z → μ + μ − and γ+jets events, as well as the p T balance between jets in dijet and multijet events [43].The uncertainty in the jet energy scale is a few percent and depends on p T and η.The impact of jet energy scale uncertainties in event yields is evaluated by varying the jet energy corrections within their uncertainties and propagating the effect to the final result.The impact of the scale uncertainties in the category yields is largest for those targeting VBF, and can be as high as 15%, but is less than 3% for the resolution.The effect of the scale uncertainties is observed in Fig. 1 (right), where the combined uncertainty has large contributions from the jet energy scale uncertainty.
• Trigger efficiency: the uncertainty in the efficiency of the trigger selection is measured with Z → e + e − events using the T&P technique.The size of its uncertainty is less than 1%.An additional uncertainty is introduced to account for a gradual shift in the timing of the inputs of the ECAL L1 trigger in the region at |η| > 2.0, which caused a specific trigger inefficiency during 2016-2017 data taking [20].
The total impact on the B(H → e + e − ) measurement from experimental systematic uncertainties is ±0.2 × 10 −4 , with the largest contributions arising from uncertainties affecting the electron energy scale.

Results
This section describes the construction of the statistical model used in the maximum likelihood fit, and presents the results of the analysis.

Signal models
The shape of the m ee distribution for simulated signal is parameterized independently for each production process entering each event category, and for each year separately.Each signal model is constructed from a sum of up to five independent Gaussian functions.To allow for the extraction of limits on B(H → e + e − ) as a function of m H , signal models are constructed from a simultaneous fit to samples at m H = 120, 125, and 130 GeV.Each parameter for each Gaussian is a linear function of m H .The optimal number of Gaussian functions is chosen by performing an F-test [49], to avoid overfitting to statistical fluctuations resulting from the limited size of the simulated samples.The sum of models for each production process is shown for the highest S/B analysis categories targeting VBF and ggH production in Fig. 3, for a nominal m H of 125 GeV.

Background models
The discrete profiling method [29] is used to account for the uncertainty related to mismodelling of the background.This method treats the choice of the function used to model the background, which is unknown a priori, as a discrete nuisance parameter.For each category, a range of analytical functions is used as candidates to fit the m ee distribution.These include generic smoothly-falling functions such as sums of exponential functions, Laurent series', and Bernstein polynomials, as well as physicsinspired models chosen to fit the dominant DY component of the background, which comprise modified Breit-Wigner functions.These functions have been used successfully to fit similar background processes in previous CMS analyses [10,28].The final set of functions considered is chosen by performing an F-test for each candidate, with a loose requirement on the goodness-of-fit.
Each candidate function is considered in the final signal-plusbackground (S+B) fit to data.The function that results in the best overall fit is chosen for each value of the parameter of interest.In the fit, the normalization and shape parameters for the background functions are allowed to vary freely.The resulting confidence intervals therefore account for possible changes in the preferred function choice.The bias on B(H → e + e − ) resulting from a particular choice of background function was also checked for each category and found to be negligible.
Finally, background events resulting from both SM H → γγ production, where the photons are misreconstructed as electrons, and from the Dalitz decay [50] (H → e + e − γ), where the photon is included in the electron reconstruction, were studied with simulated events.Each process was found to contribute around 0.1% and 0.2% to the inclusive analysis categories, respectively, for a B(H → e + e − ) at the expected limit.These processes are, therefore, neglected in this analysis.

Limits on B(H → e + e − )
A simultaneous maximum likelihood fit is performed across all analysis categories.The (S+B) m ee distributions are shown for the highest S/B analysis categories targeting ggH and VBF production in Fig. 4. Upper limits are set on B(H → e + e − ) using the CL s modified frequentist criterion [51][52][53] for m H hypotheses between 120-130 GeV.This construction uses the profile-likelihood ratio as the test statistic [54], under the asymptotic approximation.Limits for a range of m H hypotheses, along with the associated uncertainty intervals, are shown in Fig. 5.At the current best measured value of m H = 125.38GeV [55], the observed (expected) 95% CL limit is B(H → e + e − ) < 3.0 × 10 −4 (3.0 × 10 −4 ).Assuming that the only deviation from the SM prediction is from a modified Yukawa coupling to the electron, the limit translates into an upper bound on the Higgs boson effective coupling modifier to electrons of |κ e | < 240.A breakdown of the expected and observed limits on  B(H → e + e − ) is shown per analysis category in Fig. 6.Tabulated results are provided in the HEPData record for this analysis [56].

Summary
A search for the Higgs boson decaying to an e + e − pair is performed using proton-proton collision data collected at √ s = 13 TeV with the CMS experiment at the LHC between 2016-2018, corresponding to an integrated luminosity of 138 fb −1 .The analysis uses categories targeting Higgs boson production via gluon fusion and vector boson fusion, with dedicated boosted decision tree classifiers trained for each production mode to enhance the sensitivity of the resulting categories.A maximum likelihood fit to the dielectron mass distribution is performed simultaneously in each analysis category to extract an upper limit on the Higgs boson to electron pair branching fraction; the resulting observed (expected) limit at the 95% confidence level on the branching fraction for H → e + e − decays is 3.0 × 10 −4 (3.0 × 10 −4 ).This is the most stringent limit on the Higgs boson branching fraction to an e + e − pair to date.When compared with the previous best limit from the CMS Collaboration [18], where analysis categories were constructed us- ing a selection on electron and jet kinematics, the improvement in sensitivity presented in this Letter is attributed primarily to the use of BDT classifiers, which significantly improve the S/B ratio of analysis categories.Accounting both for the increase in integrated luminosity and centre-of-mass energy, the use of BDT classifiers further improves the limit on B(H → e + e − ) by a factor of approximately 1.5.

Fig. 2 .
Fig.2.Distribution of the output score of the ggH BDT in simulated background and signal events, and data (left).The ggH and VBF signals are scaled such that they are visible.Category boundaries targeting ggH Higgs boson production are denoted with dashed lines.Events with scores in the grey shaded region are discarded from the analysis.The right plot shows the distribution of the output score of the ggH BDT in a control region around the Z boson mass.Agreement is compared between the DY simulation (filled histogram) and data (black points).The combined impact of the statistical and systematic uncertainties in simulation is shown by the red shaded band, where the sources contributing to the systematic component are identical to those included in Fig.1.Residual differences between data and simulation are smaller than the theory uncertainties in the ggH cross section and event kinematics, which are included in the final maximum likelihood fit.

Fig. 3 .
Fig. 3. Signal models for the highest S/B categories targeting ggH and VBF processes, integrated over production processes, for Higgs boson events simulated at m H = 125 GeV.Contributions from each of the three years are shown by the dashed lines.The models are normalized to unit area.The σ eff is the smallest interval containing 68.3% of the m ee signal distribution.

Fig. 4 .
Fig. 4. The signal-plus-background model fit to the m ee distribution for the highest S/B analysis categories targeting the ggH (left) and VBF (right) processes.The signal model for each category is also shown, scaled to the observed limit at m H = 125.38GeV.The one (green) and two (yellow) standard deviation bands show the uncertainties in the background component of the fit.The lower panel shows the residuals after subtraction of this background component.The background functions describe the data well, with no excess observed.

Fig. 6 .
Fig. 6.Expected and observed limits on B(H → e + e − ) for each analysis category, and all categories combined The results here are computed for m H = 125.38GeV.