The Higgs boson (H) was discovered by the ATLAS and CMS experiments at the LHC in 2012 with mass m H ≈ 125 GeV [1][2][3]. Measurements of the properties of the Higgs boson, including the H decay branching fractions, are thus far found to be consistent with the expectations of the standard model (SM) [4][5][6][7][8][9][10][11]. Previous studies based on the combined results from the ATLAS and CMS experiments constrain the inclusive branching fraction of potential beyond-the-SM (BSM) undetected visible decays of the Higgs boson to be <0.12 and <0. 16 at the 95% confidence level (CL), respectively [4,5].
The lepton-flavor violating (LFV) decays H → eμ, H → eτ, or H → μτ are forbidden in the SM but may arise in BSM theories with more than one Higgs boson doublet [12,13], models with flavor symmetries [14], the Randall-Sundrum model [15][16][17][18][19], composite Higgs models [20,21], certain supersymmetric models [22][23][24], and others [25][26][27][28][29]. In these models, the LFV decays can occur through the off-diagonal LFV Yukawa couplings Y eμ , Y eτ , or Y μτ , which couple the Higgs boson with leptons of different flavor. The presence of these off-diagonal LFV Yukawa couplings may enhance processes such as μ → 3e, μ → e conversion, and μ → eγ that could proceed via a virtual Higgs boson exchange [30,31]. In particular, the most stringent limit on BðH → eμÞ is obtained indirectly from the limit on μ → eγ [32] to be <10 -8 [33]. However, the indirect limit on H → eμ assumes the SM values for the not yet tightly constrained Yukawa couplings Y μμ [34,35] and the unmeasured Y ee . For example, should Y μμ be smaller than the SM prediction, the indirect constraints on BðH → eμÞ mentioned would be loosened. It also assumes the flavor changing neutral current is dominated by the Higgs boson contribution. Hence, a direct search for H → eμ remains important. The most stringent direct limit on BðH → eμÞ, thus far, was set by the ATLAS experiment at an observed (expected) limit of 6.2ð5.9Þ × 10 -5 at 95% CL with a proton-proton (pp) collision dataset at a center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of 139 fb -1 [36].
LFV could also arise in decays of additional Higgs bosons in the Type-III two Higgs doublet model (2HDM) [13]. Recent studies have shown that searching for additional Higgs bosons with a mass below twice the W boson mass in the LFV decay channels is particularly important to constrain the Type-III 2HDM parameter space [37]. An additional Higgs boson with a mass larger than twice the W boson mass is expected to decay primarily into a W þ W - pair which dilutes the rate of LFV decays.
This paper reports a search for a LFV decay in the eμ channel of both H and of an additional Higgs boson (X) with a mass, m X , in the range 110-160 GeV. The upper range of 160 GeV corresponds to twice the W boson mass. The search is performed with data recorded by the CMS experiment in pp collisions at a center-of-mass energy of 13 TeV during the period from 2016 to 2018 and corresponding to an integrated luminosity of 138 fb -1 . The analysis in this paper is optimized for the two dominant production modes of the Higgs boson at the LHC: gluon fusion (ggH) and vector boson fusion (VBF). The final state of interest in both production modes consists of a prompt, oppositely-charged electron-muon pair. Subdominant production modes of the Higgs boson in association with a vector boson (W or Z) are not considered due to the smaller cross sections.
This paper is organized as follows: a description of the CMS detector is given in Sec. II, the collision data and simulated samples are discussed in Section III, the event reconstruction is described in Sec. IV, and the event selection is described in Sec. V. The event categorization is described in Sec. VI. Signal and background modeling, and systematic uncertainties are described in Secs. VII and VIII, respectively. Results are presented in Sec. IX, and the summary is given in Sec. X.
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. Forward calorimeters extend the pseudorapidity, η, coverage provided by the barrel and endcap detectors. Muons are detected in gasionization chambers embedded in the steel flux-return yoke outside the solenoid. 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. [38].
Events of interest are selected using a two-tiered trigger system. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of approximately 100 kHz within a fixed latency of approximately 4 μs [39]. The second level, the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing that reduces the event rate to approximately 1 kHz before data storage [40].
This search is carried out using pp collision data collected by the CMS experiment from 2016-2018 at a centerof-mass energy of 13 TeV with the integrated luminosity being 36.3 in 2016, 41.5 in 2017, and 59.8 fb -1 in 2018, respectively. Single-electron or -muon triggers with isolation criteria are used to collect the data. The transverse momentum, p T , thresholds for the electron (muon) trigger are 27 (24), 32 (27), and 32 (24) GeV in the 2016, 2017 and 2018 datasets, respectively.
Simulations are used to model the signal and background events. To model the parton showering, hadronization, and underlying event properties, PYTHIA [41] version 8.240, with the CP5 underlying event tune [42] is used in all cases. The NNPDF3.1 parton distribution functions (PDFs) are used in the simulations [43]. The simulation of interactions in the CMS detector is based on Geant4 [44]; the same reconstruction algorithms are used as for data.
The Higgs bosons are produced at the LHC predominantly via the ggH mode [45], the VBF mode [46], and in association with a vector boson (W or Z) [47]. Signal samples of H → eμ and X → eμ with a hypothesized m X of 110, 120, 130, 140, 150, and 160 GeV are generated for the ggH and VBF modes at next-to-leading order (NLO) accuracy in perturbative quantum chromodynamics (QCD) with the POWHEG v2.0 generator [48][49][50][51][52][53] using the implementation described in Refs. [54,55], interfaced with PYTHIA. The simulated Xs are assumed to have narrow width. The Herwig 7.2 generator [56] with the CH3 underlying event tune [57] interfaced with the POWHEG v2.0 generator, is used to produce alternative samples for the VBF signal. These samples are used to evaluate the systematic uncertainty in the kinematic distributions of the final state particle in VBF production due to different choices of parton shower simulation [58].
Background events from H decaying to a pair of τ leptons are simulated for all three dominant production modes at the LHC at NLO with the same POWHEG v2.0 generator as the signals, interfaced with PYTHIA. Background events from H decaying to a pair of W bosons are generated similarly for the ggH and VBF modes only as the contribution of other production modes is negligible.
The MadGraph5_aMC@NLO generator [59] (version 2.6.5) is used to simulate the single W=Z backgrounds produced by VBF in association with two or more jets from electroweak vertices (VBF W=Z þ jets) at leading order with the MLM jet matching and merging schemes [60]. Drell-Yan (DY), single W with jets from QCD vertices (QCD W þ jets), and diboson (WW, WZ, ZZ) events are simulated with the same generator at NLO, with the FxFx jet-matching and merging scheme [61]. Top quark-antiquark pair and single top quark production are generated at NLO with POWHEG v2.0.
All samples include the effects of additional pp interactions in the same or adjacent bunch crossings, referred to as pileup. The distribution of the number of pileup interactions in simulation is also weighted to match the one observed in data.
The particle flow (PF) algorithm [62] reconstructs and identifies particles in an event through an optimized combination of information from the various subdetectors of the CMS detector. The identification of the particle type (photons, electrons, muons, charged and neutral hadrons) plays an important role in determining the direction and energy of each reconstructed particle (PF candidates). The primary vertex (PV) is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone, as described in Sec. 9.4.1 of Ref. [63].
An electron is identified as a track from the PV combined with one or more ECAL energy clusters. These clusters correspond to the electron and possible bremsstrahlung photons emitted when passing through the tracker. Electrons are accepted in the range of jηj < 2.5, except for 1.44 < jηj < 1.57, the transition region between the barrel and endcap calorimeters, because the reconstruction of an electron object in this region is not optimal. Electrons with p T > 10 GeV are identified with an efficiency of 80% using a multivariate discriminant that combines observables sensitive to the amount of bremsstrahlung energy deposited along the electron trajectory, the geometric and momentum matching between the electron trajectory and the associated clusters, and the distribution of the shower energy in the calorimeters [64]. Electrons identified as originating from photon conversions are removed. The electron momentum is estimated by combining the energy measurement in the ECAL with the momentum measurement in the tracker. The momentum resolution for electrons with p T ≈ 45 GeV from Z → ee decays ranges from 1.6 to 5.0%. It is generally better in the barrel region than in the endcaps, and also depends on the bremsstrahlung energy emitted by the electron as it traverses the material in front of the ECAL [64,65].
Muons are detected in the region of jηj < 2.4 with drift tubes, cathode strip chambers, and resistive-plate chambers. Matching muons to tracks measured in the silicon tracker results in a p T resolution of 1% in the barrel and 3% in the endcaps for muons with p T up to 100 GeV. Overall, the efficiency to reconstruct and identify muons is greater than 96% [66].
The electron (muon) isolation is determined relative to its p l T values, where l is e (μ), by summing over the scalar p T of the PF particles within a cone of ΔR ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðΔηÞ 2 þ ðΔϕÞ 2 p ¼ 0.3 (0.4) around the lepton (where ϕ is azimuthal angle in radians), divided by p l T :
where p PV charged T , p neutral T , and p γ T are the p T of charged hadrons, neutral hadrons, and photons within the cone, respectively. The neutral particle contribution to the isolation from pileup, p PU T ðlÞ, is estimated for the electron from the area of jets and their median energy density in each event [67]. For the muon, half of the p T sum of the charged hadrons not coming from the PV within the isolation cone is used instead. The factor of 0.5 is estimated from simulations to be the ratio of neutral particle to charged hadron production in inelastic pp collisions [66]. The charged-particle contribution to the isolation from pileup is rejected by requiring all tracks to originate from the PV. An isolation requirement of I e rel < 0.10 (I μ rel < 0.15) is imposed to suppress backgrounds of jets misidentified as an electron (muon).
Charged hadrons are identified as charged particle tracks neither identified as electrons, nor as muons. Neutral hadrons are identified from HCAL energy clusters not assigned to any charged hadron or from an excess in the combined ECAL and HCAL energy with respect to the expected charged-hadron energy deposit.
For each event, hadronic jets are clustered from these reconstructed particles using the infrared and collinear safe anti-k T algorithm [68,69] with a distance parameter of 0.4. Jet momentum 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. Pileup can contribute additional tracks and calorimetric energy depositions to the jet momentum. To mitigate this effect, charged particles identified to be originating from pileup vertices are discarded and an offset correction is applied to correct for remaining contributions. Jet energy corrections are derived from simulation to bring the measured response of jets to that of particle level jets on average. In situ measurements of the momentum balance in dijet, photon þ jet, Z þ jet, and multijet events are used to account for any residual differences in the jet energy scale between data and simulation [70]. The jet energy resolution amounts typically to 15%-20% at 30 GeV, 10% at 100 GeV, and 5% at 1 TeV [70]. Additional selection criteria are applied to each jet to remove jets potentially dominated by anomalous contributions from various subdetector components or reconstruction failures. Jets are required to have a p T > 30 GeV, jηj < 4.7, and be separated from each lepton of the identified eμ pair by ΔR > 0.4. Jets originating from b hadron decays and detected within the tracker acceptance of jηj < 2.5 are tagged using a deep neural network based algorithm, DeepJet, using a working point with a 94% b-jet identification efficiency at a 10% misidentification rate for lightflavor quark and gluon jets in t t events [71].
Hadronic τ decays (τ h ) are reconstructed from jets, using the hadrons-plus-strips algorithm [72], which combines 1 or 3 tracks with energy deposits in the calorimeters, to identify the tau decay modes. To distinguish genuine τ h decays from jets originating from the hadronization of quarks or gluons, and from electrons, or muons, the DeepTau algorithm is used [73]. Information from all individual reconstructed particles near the τ h axis is combined with properties of the τ h candidate and the event. The rate at which jets are misidentified as τ h by the DeepTau algorithm depends on the p T and whether it was initiated by a quark or gluon. A working point with an 80% τ h identification efficiency and a 0.05-0.95% misidentification rate for jets is used.
The missing transverse momentum, ⃗p miss T , is computed as the negative vectorial p T sum of all the PF candidates in an event [74], with its magnitude labeled as p miss T . Corrections to the reconstructed jet energy scale are propagated to the ⃗p miss T . Anomalous high-p miss T events can originate from various reconstruction failures, detector malfunctions, or noncollision backgrounds. Such events are rejected using event filters designed to identify more than 85%-90% of the spurious high-p miss T events with a misidentification rate of less than 0.1% [74].
The signal topology consists of an oppositely charged electron-muon pair with possible additional jets. Events with an oppositely charged electron-muon pair separated by ΔR > 0.3 are selected. Both the electron and muon are required to have a longitudinal and a transverse impact parameter within 5 and 2 mm from the PV, respectively.
The invariant mass of the eμ pair, m eμ , is required to fall in the range of 100-170 GeV such that signals with the lowest (highest) m X ¼ 110ð160Þ GeV targeted in this search are fully contained. The m eμ window is intentionally chosen to lie beyond the peak of the t t background distribution, thus selecting a region where it falls smoothly. Backgrounds from H → ττ and H → WW also peak below the mass window since part of the H four-momentum is carried away by the final-state neutrinos.
The p T of the electron (muon), p e T (p μ T ), collected by the single-electron (single-muon) triggers is required to be larger than 29 (26) GeV in 2016, 34 (29) GeV in 2017, and 34 (26) GeV in 2018. These p T requirements are chosen to be slightly above the p T thresholds of the triggers so that the efficiency of the triggers is nearly 100%. For electrons (muons) that do not pass the single-electron (single-muon) trigger requirements, their p T are required to be larger than 25 (20) GeV in all years. Events containing additional reconstructed electrons, muons, or hadronically decaying tau candidates are vetoed. Events with at least one b-tagged jet are also vetoed to suppress the t t and single top quark backgrounds.
Events are first divided into two broad categories to enhance the signal from either the ggH or the VBF production mechanisms. Events with two or more jets where the two highest p T jets have an invariant mass m j 1 j 2 > 400 GeV and a pseudorapidity separation jΔηðj 1 ; j 2 Þj > 2.5 are classified as the VBF production category. Otherwise, events enter the ggH production category. The m eμ distributions of the data, the simulated backgrounds, and signals of H → eμ are shown in Fig. 1 for both categories. The QCD multijet background shown is estimated from a control region of the data using events with an eμ pair of the same electric charge and extrapolated to the signal region as a function of jet multiplicity and the ΔR separation of the eμ pair as described in Ref. [75]. The data and background simulations show good agreement within the statistical and the SM cross section uncertainties combined.
The two broad categories are further split according to the signal purity using the output of boosted decision trees FIG. 1. The m eμ distributions of the data, simulated backgrounds and signals of H → eμ in the ggH (upper) and the VBF categories (lower). A BðH → eμÞ ¼ 0.2% is assumed for the signal for illustration. The lower panel in each plot shows the ratio of the data to the total estimated background. The uncertainty band corresponds to the background uncertainties, adding in quadrature the statistical and the SM cross section uncertainties.
(BDTs) trained with the XGBoost package [76]. The BDTs are trained separately for the ggH and the VBF categories. The BDT discriminants range from 0 for backgroundlike events to 1 for signal-like events. For both BDTs, a mixture of simulated signal events is used in the training including events of H → eμ and X → eμ at m X ¼ 110, 120, 130, 140, 150, and 160 GeV from both the ggH and the VBF production modes. Kinematics variables from the dominant sources of backgrounds of di-leptonic decays of t t and WW diboson events are used in the training. All events used in the training are from Monte Carlo (MC) simulations described in Sec. III. The simulated signals of H → eμ and backgrounds are weighted according to their expected yields from the SM cross sections. The simulated signal samples of X → eμ are weighted according to their relative SM-like production cross sections as evaluated in Ref. [77] as a function of m X . Their total weights are matched to that of the backgrounds in the training to ensure the larger total weights of the background samples does not lead to BDTs with poor signal identification efficiency. Each signal event is additionally reweighted by the inverse of its expected mass resolution during training. The mass resolution is the uncertainties of m eμ propagated from the expected uncertainty of the lepton p T measurements. This reweighting allows the BDTs to assign more importance in classifying signal events with high mass resolution. The ggH and VBF BDT discriminant distributions of the data, the simulated backgrounds, and signals of H → eμ are shown in Fig. 2 for the ggH and VBF category, respectively. The data and background simulations show good agreement within the statistical and the SM cross section uncertainties combined.
The BDT input variables are chosen such that the BDTs do not make use of m eμ to discriminate between the signal and background. This ensures background events with m eμ close to the signal resonance are not preferrentially assigned a higher BDT discriminant which distorts their smoothly falling shape to form spurious signal resonances along m eμ . For example, the eμ system's p T scaled by m eμ , p eμ T =m eμ , is used instead of p eμ T which is correlated with m eμ . The background samples are also reweighted to match the shape of the m eμ distribution of the signals during training to further ensure that the BDTs do not benefit from using m eμ to discriminate between the signal and background. No sculpting of the m eμ distribution is observed for the MC background events in different ranges of the BDT discriminants. The BDT discriminant distributions of the simulated signals at different Higgs boson masses are also observed to be similar. The distributions of p miss T , which is the most discriminating variable in both the ggH and VBF category, are shown in Fig. 3 for both categories.
The input variables to the ggH BDT include the absolute pseudorapidities of the electron, jη e j, and of the muon, jη μ j, the ratio of the eμ system's p T to its invariant mass, p eμ T =m eμ , and the pseudorapidity separation of the eμ pair, A BðH → eμÞ ¼ 1.0% is assumed for the signal for illustration. The lower panel in each plot shows the ratio of the data to the total estimated background. The uncertainty band corresponds to the background uncertainties, adding in quadrature the statistical and the SM cross section uncertainties. Vertical lines in the plots illustrate boundaries of the subcategories: ggH cat 0-3 and VBF cat 0-1, as defined in Sec. VI B. Events in the shaded region of the VBF category with a VBF BDT discriminant less than 0.78 are discarded since their sensitivity is an order of magnitude lower than other subcategories.
jΔηðe; μÞj. Variables related to the ⃗p miss T including p miss T and its azimuthal separation to the eμ system, Δϕð ⃗p miss T ; p eμ T Þ, are also used to discriminate the neutrinoless LFV decay against backgrounds with neutrinos in the final state. The number of jets in each event is also added as an input variable.
Additional jet variables are added for events with at least one jet, including the p T of the leading jet, p j 1 T , the absolute pseudorapidity of the leading jet, jη j 1 j, and the pseudorapidity separation of the leading jet to the eμ system, jΔηðj 1 ; eμÞj. For events with at least two jets, the scalar p T of all jets is added. Observables sensitive to the angular and p T correlations between the eμ system and the two highest p T jets are also included, including the p T -balance ratio:
and the p T -centrality:
If jets are absent in an event, the undefined jet variables are handled by the sparsity-aware split finding algorithm in the XGBoost package [76], with the exception of p j 1 T set to be zero in events with no jets. When jet variables are used at a decision split of a tree, the sparsity-aware algorithm assigns events with an undefined value to the direction that minimizes the loss function.
The input variables to the VBF BDT are the same to that of the ggH BDT with a few exceptions: Δϕð ⃗p miss T ; p eμ T Þ, jη j 1 j, and jΔηðj 1 ; eμÞj are dropped due to their insignificant contributions to the VBF BDT training. Instead, the Zeppenfeld variable [78], defined as
is added along with m j 1 j 2 and jΔηðj 1 ; j 2 Þj.
The ggH and VBF categories are further split according to the ggH and VBF BDT discriminant value to optimize the expected sensitivity of the search. The expected sensitivity is estimated from the expected significance of discovery in the asymptotic approximation [79] from a signal-plus-background (S þ B) fit to the m eμ distribution in the data within 100-170 GeV, overlayed with a simulated signal of BðH → eμÞ ¼ 5.9 × 10 -5 , the most stringent direct expected limit of BðH → eμÞ up-to-date [36]. In these fits, the signal peaks are modeled with a sum of three (two) Gaussians for the ggH (VBF) production signals in the ggH categories. For both the ggH and VBF production signals in the VBF categories, a sum of two Gaussians is used. The number of Gaussians chosen are motivated by the likelihood-ratio test [80], as explained in Sec. VII. The total expected background is modeled from data directly with a third (first) order Bernstein polynomial for the ggH (VBF) category. Subcategory boundaries are determined separately in the ggH and the VBF categories by iteratively scanning in steps of 0.01 for a cutoff along the ggH and the VBF BDT discriminants, respectively, that maximizes the FIG. 3. The p miss T distributions of the data, simulated backgrounds and signals of H → eμ in the ggH (upper) and the VBF categories (lower). A BðH → eμÞ ¼ 1.0% is assumed for the signal for illustration. The lower panel in each plot shows the ratio of the data to the total estimated background. The uncertainty band corresponds to the background uncertainties, adding in quadrature the statistical and the SM cross section uncertainties. total expected sensitivity. This procedure is repeated until the further gain in sensitivity is less than 1%.
Four optimized subcategories are defined for the ggH category, named as "ggH cat 0," "ggH cat 1," "ggH cat 2," and "ggH cat 3," which correspond to events of decreasing ggH BDT discriminant of 0.89-1.00, 0.77-0.89, 0.46-0.77, and 0.00-0.46, respectively. Similarly, three optimized subcategories are defined for the VBF category, "VBF cat 0," "VBF cat 1," and "VBF cat 2," corresponding to events with a VBF BDT discriminant between 0.94-1.00, 0.78-0.94, and 0.00-0.78 respectively. Events from the least sensitive category "VBF cat 2" are discarded. Table I summarizes the definition, the expected background (B), and signal yield of H → eμ at B ¼ 10 -4 (S) in each categories at an integrated luminosity of 138 fb -1 . An estimate of the expected significance in each category by S= ffiffiffi ffi B p is also listed. The yields are estimated by the number of MC events within a m eμ interval of 125 GeV AE σ eff , where σ eff is half of the smallest symmetric interval that contains 68% of the signal events in each category.
The m eμ distributions of simulated signal events are fit with a sum of Gaussian distributions for each production mode, category, and mass of the Higgs boson. The number of Gaussians is chosen with the likelihood ratio test [80], such that the next higher order does not give a significantly better fit at a p-value of 0.05. A sum of three Gaussians is determined to be sufficient for the signals from the ggH production mode in the ggH category, while a sum of two Gaussians is sufficient for the rest. When carrying out the fits, the means are fit as a sum of the known simulated m X or m H and a small floating shift due to initial/final-state radiations and detector effects. Example fits of the signal models to the simulated H → eμ signal are shown in Fig. 4 for the analysis categories ggH cat 0 and ggH cat 3, as well as VBF cat 0 and VBF cat 1, summing events from both the ggH and VBF production modes. σ eff for each distribution is included as an illustration of the signal resolution. Example fits of the signal models to the simulated H → eμ signal in the analysis categories ggH cat 0 and ggH cat 3 (upper), as well as VBF cat 0 and VBF cat 1 (lower), summing events from both the ggH and VBF production modes. Half of the smallest symmetric m eμ interval that contains 68% of the signal events, σ eff , is shown in the legends for each signal as an illustration of the signal resolution. The signal resolution in general improves with the signal purity of the analysis categories. The signal resolution in general improves with the signal purity of the analysis categories since signal events are reweighted by the inverse of their mass resolution during training of the BDTs as mentioned in Sec. VI. The m eμ distributions of a Higgs boson with mass between the simulated mass points are interpolated by fitting the parameters and normalizations of the sum of Gaussians with second-order polynomials as a function of the Higgs boson mass. The background in each category is modeled with a Bernstein polynomial. Orders of the polynomials are chosen with a bias study as follows. The m eμ distribution in data from 100-170 GeV is first fit with three distinct functional forms: a Bernstein polynomial, a sum of exponential functions, and a sum of power law functions. An optimal order for each function is chosen with the likelihood-ratio test [80] at a p-value of 0.05. Then, ensembles of 2000 pseudo-experiments are generated with the m eμ distributions drawn from each of the three background models, with or without an injection of a signal at the simulated m X points with a branching fraction of 10 -4 . The pseudoexperiment mass spectra are fit with a Bernstein polynomial with an order equal to or higher than the chosen order in the first step. The signal yield from these fits would in general differ from the injected yield since different background models are used to generate and fit the signal peaks. The bias of a model choice is evaluated as the average difference of the fit signal yield to the generated yield divided by the uncertainty in the fit yield in the pseudoexperiments. The final order of the Bernstein polynomial used to model the background in each category is then chosen by requiring the bias to be less than 20% across all generating functions and ensembles of pseudoexperiments. The third order is chosen for all ggH subcategories, while the second and the first order are chosen for "VBF cat 0" and "VBF cat 1," respectively.
The systematic uncertainty associated to the bias of the background model choice is modeled by adding a signallike background shape to the background models. The signal-like background shape is drawn directly from the signal models in each category. The normalization of the signal-like background is implemented as a nuisance parameter modeled with a Gaussian constraint of zero mean and a standard deviation equal to the maximum of the pseudo-experiment averaged signal yield fit over the three background models in the bias study with no signal injected, as described in Section VII. The maxima are no more than 20% of the statistical uncertainties in the fits. The standard deviations amount to a Higgs-like signal yield with a branching fraction BðH=X → eμÞ of 0.4-2.9 × 10 -5 across the categories. This is a dominant source of systematic uncertainty, contributing 6.9-14.4% of the total uncertainty on the best fit of the signal yield, depending on m H =m X . Besides the systematic uncertainty associated to the bias of the background model and the statistical uncertainty of the fits, there are no additional systematic uncertainties in the background models as they are derived directly from data.
The simulated signals are affected by various sources of experimental and theoretical systematic uncertainties. These uncertainties affect both the yield and the shape of the m eμ distributions. The systematic uncertainties are incorporated as nuisance parameters in the S þ B likelihood fit of the m eμ distribution. Log-normal constraints are assumed for uncertainties affecting the yield, and Gaussian constraints are assumed for uncertainties affecting the fit parameters of the m eμ distribution. The uncertainties affecting the yield have negligible effects on the signal shapes in general. All the uncertainties are treated as correlated between the categories, except for systematic uncertainties from the interpolation of signal shapes. The list of yield uncertainties is summarized in Table II for the ggH and VBF production modes separately.
The uncertainties in the electron (muon) momentum scale and resolution affect the means and widths of the signal models. These uncertainties are measured in Z → ee (Z → μμ) events in data and simulation in the H → ZZ → 4lðl ¼ e; μÞ analysis [81]. They are estimated to be 0.1% for the means and 10.0% for the widths of the signal models.
The uncertainties in the reconstruction, single-lepton trigger, offline identification, and isolation efficiencies of electrons and muons are respectively measured in Z → ee and Z → μμ events with the "tag-and-probe" method [82] in data and simulated events. They amount to be 1.8-2.6% for electrons and 0.2-0.4% for muons [64,66]. The lepton identification and isolation uncertainties are treated as correlated between the data-taking years, while the trigger uncertainties are treated as uncorrelated.
The uncertainties in the jet energy scale and resolution from different sources are evaluated as functions of the jet p T and η [70]. Jets with p T < 10 GeV are classified as unclustered energy. The uncertainties in the unclustered energy scale for charged particles, neutral hadrons, photons, and very forward particles are evaluated separately according to the resolution of the different sub-detectors. The combined uncertainty of the unclustered energy scale is then propagated to the ⃗p miss T . Uncertainties on jets and ⃗p miss T affect both the ggH and VBF BDTs, which are used to define the categories. They are transformed into signal yield uncertainties per category, which in turn enter as nuisance parameters in the likelihood fit. The efficiency to identify a b-tagged jet with the DeepJet algorithm is different in data and simulations and affect the b-tagging veto. Scale factors dependent on the jet p T and η are applied to correct the simulation [83,84]. The uncertainties in these scale factors are taken into account.
The theoretical uncertainties in the renormalization and factorization scales, the choice of PDFs, and the value of the strong coupling constant, α S , evaluated at the Z boson mass, affect the measurement of the Higgs boson production cross sections [77]. These uncertainties in turn affect the expected signal yield and are treated as correlated between the datataking years. The QCD scales variations lead to 3.9-8.0% and 0.2-0.5% of uncertainty in the ggH and VBF cross sections, respectively, while changes in the PDFs and α S result in 3.0-3.2% and 1.9-2.1% uncertainties, respectively. The uncertainties in the event acceptance in each category due to the scales, PDFs, and α S are also taken into account. An additional uncertainty in the VBF parton shower model is assigned as the signal yield difference between the dipole shower in PYTHIA and the alternative angular-ordered shower in Herwig. This amounts to 1.9-11.4% uncertainties across the categories.
The integrated luminosities for the 2016-2018 datataking years have 1.2-2.5% individual uncertainties [85][86][87], while the overall uncertainty for the 2016-2018 period is 1.6%. They affect the overall yield of the signal expected from simulations. The uncertainty on the number of pileup vertices is evaluated by varying the pileup correction weights applied to the simulation. The variation of weight is obtained through a AE4.6% change to the total inelastic cross section at a nominal value of 69.2 mb used to estimate the pileup effect on data. The pileup uncertainties are treated as correlated between the years.
During the 2016 and 2017 data-taking periods, a gradual timing shift of the signals from the ECAL first-level trigger caused a specific trigger inefficiency in the region of jηj > 2.0. For events containing an electron with p T > 50 GeV or a jet with p T > 100 GeV, in the region of 2.5 < jηj < 3.0, the efficiency loss is between 10.0%-20.0%, dependent on p T , η, and time. Scale factors are computed to correct the detector acceptance in simulations to reflect this effect in the data. The uncertainty due to this correction is 0.1%-0.5% and is treated as correlated between the two years.
No excess of data above the background prediction has been observed for the LFV decay of H → eμ. An upper limit on the branching fraction of the decay is computed using the CL s criterion, with the profile likelihood as the test statistic, which is assumed to be distributed as in the asymptotic approximation [79,88,89]. The observed (expected) upper limit on BðH → eμÞ is 4.4ð4.7Þ × 10 -5 at 95% CL. A breakdown of the upper limit on BðH → eμÞ is shown per analysis category, and for the combination of all analysis categories is illustrated graphically in Fig. 5 and listed in Table III. Tabulated results are provided in the HEPData record [90].
The upper limit on BðH → eμÞ is also interpreted as a constraint on the LFV Yukawa couplings Y eμ [33]. The LFV decay arises at tree level from the BSM Yukawa coupling, Y eμ . The decay width ΓðH → eμÞ can be written in terms of the Yukawa coupling as,
The branching fraction BðH → eμÞ assuming H → eμ is the only BSM contribution is given by,
FIG. 6. Constraints on the lepton-flavor violating Yukawa couplings, jY eμ j and jY μe j. The observed (expected) limit in black (red) line is derived from the limit on BðH → eμÞ in this analysis. The green (yellow) band indicates the one (two) standard deviation uncertainty in the expected limit. The hashed region is excluded by this direct search. Other shaded regions represent indirect constraints derived from the null searches for μ → 3e (gray) [92], μ → e conversion (light blue) [93], and μ → eγ (dark green) [32]. The flavor-diagonal Yukawa couplings, jY ee j and jY μμ j, are assumed to be at their SM values in the calculation of these indirect limits. The purple line is the theoretical naturalness limit of jY eμ Y μe j ≤ m e m μ =v 2 , where v is the vacuum expectation value of the Higgs field. Dotted lines represent the corresponding constraints on jY eμ j and jY μe j at upper limits on BðH → eμÞ at 10 -5 ; 10 -6 ; 10 -7 , and 10 -8 , respectively.
The decay width of H is assumed to be Γ SM ¼ 4.1 MeV at m H ≈ 125 GeV [91]. The observed (expected) upper limit on the Yukawa coupling is evaluated to be ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi jY eμ j 2 þ jY μe j 2 q < 1.9ð2.0Þ × 10 -4 at 95% C.L. The result is illustrated in Fig. 6.
The observed (expected) upper limit at 95% CL on σðpp → X → eμÞ is plotted as a function of the hypothesized m X in the range 110-160 GeV on the left in Fig. 7, assuming the relative SM-like production cross sections of the ggH and VBF production modes as evaluated in Ref. [77]. An excess of events over the background-only hypothesis is observed at m X ≈ 146 GeV. The corresponding S þ B fit combining all categories is shown in Fig. 8, where events in each category are weighted by S=ðS þ BÞ, where S and B are the fit number of signal and background events in that category. The observed and expected upper limits on σðpp → Xð146Þ → eμÞ at 95% CL per-category and combined are listed in Table IV and illustrated graphically in Fig. 9. The best fit of σðpp → Xð146Þ → eμÞ combining all analysis categories is 3.89 þ1.11 -1.08 ðstatÞ þ0.57 -0.34 ðsystÞ fb, with the uncertainties dominated by the statistical uncertainties of the data. The best fit of σðpp → Xð146Þ → eμÞ per-category and combined with the corresponding local significance are also summarized in Table IV. Tabulated results are provided in the HEPData record [90]. Such an excess, however, was not reported in a search of similar sensitivity for H → eμ carried out by the ATLAS experiment which covered the m eμ range of the excess [36].
Searches for the lepton-flavor violation decay of H and X with a m X in the range 110-160 GeV have been performed in the eμ final state in data collected by the CMS experiment. The data correspond to an integrated luminosity of 138 fb -1 of pp collisions at a center-of-mass energy of 13 TeV. The observed (expected) upper limit on the branching fraction of the H decay BðH → eμÞ is found to be 4.4ð4.7Þ × 10 -5 at 95% confidence level, which is the most stringent direct limit set thus far. Upper limits on the cross sections of pp → X → eμ are set in the m X range 110-160 GeV at 95% confidence level. This is the first result of a direct search for X → eμ, with m X below twice the W boson mass. The largest excess of events over the expected background is observed with a local (global) significance of 3.8 (2.8) standard deviations at an invariant mass of the eμ final state of around 146 GeV. The observed significance of this excess is insufficient to draw any conclusions. More data will be needed to clarify the nature of the excess.