Search for nonresonant pair production of highly energetic Higgs bosons decaying to bottom quarks

A search for nonresonant Higgs boson (H) pair production via gluon and vector boson (V) fusion is performed in the four-bottom-quark final state, using proton-proton collision data at 13 TeV corresponding to 138 fb$^{-1}$ collected by the CMS experiment at the LHC. The analysis targets Lorentz-boosted H pairs identified using a graph neural network. It constrains the strengths relative to the standard model of the H self-coupling and the quartic VVHH couplings, $\kappa_{2V}$, excluding $\kappa_{2V}$ = 0 for the first time, with a significance of 6.3 standard deviations when other H couplings are fixed to their standard model values.


1
The discovery of a Higgs boson by the ATLAS and CMS Collaborations at the CERN LHC in 2012 [1][2][3] started a new era in experimental high energy physics. Since its original observation, significant efforts have been made to measure the properties of the Higgs boson, including its self-coupling. The production of a pair of Higgs bosons (HH) is a rare process that provides access to direct measurements of the Higgs boson trilinear self-coupling (λ) and the quartic couplings between two Higgs bosons and two vector bosons.
The HH production takes place primarily via gluon-gluon fusion (ggF) and vector boson fusion (VBF). In the standard model (SM), the ggF cross section at √ s = 13 TeV for a Higgs boson mass of 125 GeV, calculated at next-to-next-to-leading order (NNLO) precision in perturbative quantum chromodynamics (QCD), is 31.1 +2. 1 −7.2 fb [4,5] and is sensitive to the values of λ and the top quark Yukawa coupling, y t . Variations from their SM values are parametrized as κ λ = λ/λ SM and κ t = y t /y SM t . The latter has been measured to be consistent with the SM [6, 7]. The VBF cross section, calculated at next-to-NNLO (N 3 LO) precision, is 1.726 ± 0.036 fb [8] and depends on the strength of the self-coupling and on the interaction of a pair of W or Z bosons (jointly denoted as V) with a single Higgs boson (VVH) or a pair of Higgs bosons (VVHH), whose values with respect to the SM predictions are parametrized as κ V and κ 2V , respectively. The ATLAS and CMS Collaborations have performed studies of HH production at √ s = 7, 8, and 13 TeV in the separate bb γγ [9-12], bb τ − τ + [13,14], bbbb [15][16][17][18][19][20], and bbVV [21-24] channels, as well as in combinations of channels [25][26][27].
This Letter reports the results of a search for nonresonant HH production via both the ggF and VBF modes in the channel where each Higgs boson decays to a bottom quark-antiquark pair (bb). It is based on data from proton-proton collisions at the LHC at √ s = 13 TeV, collected by the CMS experiment in 2016-2018, corresponding to an integrated luminosity of 138 fb −1 . We select events with both Higgs bosons in the highly Lorentz-boosted regime, i.e., with sufficiently large transverse momentum (p T ) for the decay products of each Higgs boson to become merged into a single large-radius jet. Events are separated into mutually exclusive categories targeting either ggF or VBF HH production, with the latter requiring a characteristic VBF signature of two additional small-radius jets in opposite forward regions. The ggF categories are designed to select SM ggF HH signal events, whereas the VBF categories target non-SM coupling hypotheses that can dramatically increase the VBF HH cross section, especially for highly Lorentz-boosted Higgs boson pairs [28][29][30][31]. Such non-SM couplings can arise for example in the models where the Higgs boson is a pseudo-Goldstone boson from a new TeVscale strong interaction [30,32], and a case of particular interest is κ 2V = 0, corresponding to vanishing VVHH couplings. The two dominant background processes, namely QCD multijet production, consisting uniquely of jets produced through the strong interaction, and top quark-antiquark pair (tt) production, are estimated utilizing dedicated data samples enriched in these processes. One of the main challenges for this search is the efficient reconstruction of the decay of a pair of Higgs bosons, while rejecting jets originating from light-flavor quarks and gluons. For this purpose, we use a novel graph neural network (GNN) algorithm, Parti-cleNet [33], to select large-radius jets resulting from H → bb decays, increasing signal purity in a mass-independent way. The search complements that in Ref. [20], which targets the same bbbb final state in the non-boosted regime and covers a region of signal phase space that is largely disjoint from that of the present analysis. Simulated samples for the ggF process are generated at next-to-leading order (NLO) precision using , and corrected as a function of HH mass based on Ref. [54]. The VBF HH samples are generated at leading order (LO) precision using MADGRAPH5 aMC@NLO 2.6.5 [55]. A range of samples corresponding to different combinations of the κ V , κ 2V , and κ λ coupling modifiers is generated. Samples for other coupling combinations are constructed as linear combinations of the original generated samples by applying appropriate event weights as described in Ref. [20]. The ggF samples are normalized to the NNLO cross sections [4,5,[56][57][58][59][60][61] corresponding to the coupling values considered. The VBF sample with SM couplings is normalized to the N 3 LO cross section [8], and the same N 3 LO/LO correction is applied to the VBF samples with modified couplings.
Samples for the tt background process are generated at NLO precision using POWHEG v2.0 [45-47, 62, 63] with the FxFx jet matching and merging scheme [64], and normalized to the theoretical cross section calculated at NNLO precision using TOP++ v2.0 [65]. The differential tt cross section as a function of the top quark p T is corrected to the NNLO QCD + NLO electroweak precision [66]. Samples for other background processes, such as the QCD multijet, single top quark, W and Z single-boson and diboson production processes, as well as single Higgs boson production in all relevant production modes, are also generated. Parton showering, hadronization, and the underlying event are modeled by PYTHIA 8.205 [67] with parameters set by the CUETP8M1 [68] and CP5 tunes [69] used for samples simulating the 2016 and 2017-2018 conditions, respectively. The NNPDF 3.0 [70] and 3.1 [71] parton distribution functions (PDFs) are used in the generation of all simulated samples. The GEANT4 [72] package is used to model the response of the CMS detector, and simulated minimum bias interactions are mixed with the hard interactions in simulated events to model additional proton-proton interactions within the same or nearby bunch crossings (pileup). The simulated events are weighted to match the pileup distributions measured in data.
The reconstructed particles are clustered into jets using the anti-k T algorithm [73] implemented in the FASTJET package [74] with a distance parameter of 0.4 (small-radius jets) or 0.8 (largeradius jets). To mitigate the effect of pileup on small-radius jets, the charged-hadron subtraction algorithm [40] is applied. For large-radius jets, the pileup-per-particle identification algorithm [75,76] is used to assign a weight to each particle prior to jet clustering based on the likelihood of the particle originating from the PV. Further corrections are applied to the jet energy as a function of the jet pseudorapidity (η) and p T to account for detector response nonlinearities [42]. The missing transverse momentum vector p miss T is computed as the negative vector sum of the transverse momenta of all the reconstructed particles in an event, and its magnitude is denoted as p miss Electrons are reconstructed and their momentum is estimated by combining the momentum measurement from the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track [37]. Muons are reconstructed as tracks in the central tracker, consistent with either a track or several hits in the muon chambers, and their momentum is obtained from the curvature of the corresponding track [38]. Electron candidates are reconstructed within the tracker acceptance of |η| < 2.5, whereas muon candidates are required to be within the muon chamber acceptance of |η| < 2.4. The electron and muon candidates are required to have p T > 10 GeV. Only tracks consistent with the PV are associated with the electrons or muons, and additional identification criteria [38,77] are applied to improve the purity of the lepton selection. Isolation criteria [78] are used to distinguish leptons originating at the PV from those originating from bottom or charm hadron decays.
Events are required to contain at least two large-radius jets, and the two with the highest p T are considered as Higgs boson candidates. A combination of several trigger algorithms is used, all requiring the total hadronic transverse energy in the event (H T ) or jet p T to be above a given threshold. After removing the remnants of soft radiation [79], a minimum triggering threshold on the jet mass is imposed in order to reduce the H T or p T thresholds relative to a trigger algorithm with no jet mass requirement but of an equivalent rate. The trigger efficiency varies between 10 and 95% for jets with p T in the range from 300 to 450 GeV and is fully efficient for jets with p T > 500 GeV.
The GNN algorithms [33,[80][81][82] have been shown to achieve state-of-the-art performance in jet classification [83], and this analysis is among the first to use the ParticleNet GNN algorithm [33] to discriminate between H → bb and QCD-induced jets. Its inputs comprise 20 measured properties of jet constituent particles, such as the p T , the angular separation between the particle and the jet axis, and the transverse and longitudinal impact parameters relative to the PV. The inputs also include eleven properties of the secondary vertices [84] within the jet cone, such as the displacement, p T , mass, and the number of associated tracks. Given the relatively large masses, long lifetimes, and high-p T and displaced decay products of b hadrons, this low-level tracking and vertexing information is important for identifying jets containing b hadron decays [85]. The algorithm treats each jet as an unordered set of its constituents, considers their correlations at different scales to extract jet features, and returns an output corresponding to the probability that the jet belongs to either the H → bb or the QCD multijet class. As shown in Ref. [86], algorithms like this one that exploit the correlations between many low-level input features outperform those that are based on fewer high-level discriminating variables for H → bb identification. To achieve independence from the jet mass, the classifier is trained using dedicated simulated samples containing spin-0 particles with a uniform mass spectrum between 15 and 250 GeV and decaying to quark-antiquark pairs, as well as a QCD multijet sample [87]. In addition, jets in each sample are reweighted to obtain uniform distributions in both p T and mass for the training. The ParticleNet algorithm yields significant improvements both in terms of performance and jet mass decorrelation compared to previous jet classifier algorithms [86,87]. The discriminant (D bb ) is calibrated using data and simulated samples dominated by QCD multijet events [88]. A boosted decision tree (BDT) classifier is used to select a sample of jets enriched in gluon splittings to bb that yields a D bb distribution similar to the HH signal. A fit to data is performed to extract p T -dependent simulation-to-data correction factors, typically ranging from 0.9 to 1.2. They are applied to correct large-radius jet selection efficiencies in the signal events. Three working points (WPs) based on D bb are used: these tight, medium, and loose WPs correspond to H → bb selection efficiencies of approximately 60, 80, and 90%, and QCD background misidentification rates of approximately 0.3, 1, and 2%, respectively.
The search for Higgs boson pairs relies on the ability to reconstruct accurately the masses of the two Higgs boson candidates. The soft-drop (SD) algorithm [89] with angular exponent β = 0 and soft radiation fraction z = 0.1, also known as the modified mass-drop algorithm [90], is usually applied to the Higgs boson jet candidate to remove soft and wide-angle radiation. The SD algorithm occasionally over-subtracts genuine jet constituents that account for a large frac-tion of the jet momentum, resulting in a jet mass (m SD ) that is negligible even for large-mass resonances. To improve the jet mass estimation, we introduce a regression algorithm based on the ParticleNet GNN architecture to predict the jet mass [91]. To train the algorithm, the target mass value that the algorithm aims to predict is chosen as the true resonance mass for the heavy spin-0 particles, or the SD mass of particle-level jets for QCD multijet events. Particle-level jets are defined by clustering stable simulated final-state particles, excluding neutrinos, using the same procedure as for reconstructed jets, and they are geometrically matched to the reconstructed jets. The training is performed with the same samples and a similar event weighting scheme as for jet classification. Data control samples, dominated by tt events, are used to calibrate the jet mass from the regression (m reg ). The small corrections to mass scale (<1%) and resolution (≈3%) are then applied to all simulated samples in the rest of the analysis, accompanied by systematic uncertainties that account for the differences in regression performance between the jets originating from W and Higgs boson decays.
Events containing at least two large-radius jets with p T > 300 GeV and |η| < 2.4, and no isolated electrons or muons, are selected and grouped into either the ggF or the VBF categories.
For the VBF categories, we require the two large-radius jets to have m reg between 50 and 200 GeV. The search region (SR) is defined as the subset of events with m reg of the p T -leading (p T -subleading) jet in the range 110-150 (90-145) GeV, while the remaining events are used as part of the background estimation method described below. As the VBF process is characterized by the presence of two forward jets with a large dijet invariant mass and a gap in η, we also require two additional small-radius jets, referred to as VBF jet candidates, with p T > 25 GeV and |η| < 4.7. They are required to be separated from both Higgs boson candidate jets by a minimal distance of 1.2 in the η-φ plane, where φ denotes the azimuthal angle. If more than two such jets are found, the two with the highest p T are selected. The absolute difference in pseudorapidity (∆η VBF jj ) between the two VBF jets must be larger than 4.0, and the invariant mass of the two VBF jets (m VBF jj ) must be larger than 500 GeV. To increase the sensitivity for the highly boosted Higgs boson pairs from VBF production with non-SM couplings, the p Tleading (p T -subleading) large-radius jet is required to have p T > 500 (400) GeV. Additionally, the azimuthal angle between the two large-radius jets (∆φ j 1 j 2 ) is required to be larger than 2.6, and the absolute difference in pseudorapidity ∆η j 1 j 2 must be less than 2.0. Three VBF event categories are defined based on the D bb scores of the two Higgs boson candidate jets. Events in the high-purity (HP) category must have both jets passing the tight D bb WP; events not in the HP category are categorized as medium purity (MP) if both jets pass the medium D bb WP; and events not in the HP or MP categories are categorized as low purity (LP) if both jets pass the loose D bb WP. The definitions of the three categories based on the corresponding D bb WPs were optimized to provide the best expected significance for the VBF HH signal hypotheses of interest such as κ 2V = 0 (with other couplings at their SM values), while enabling a robust background estimation and reliable D bb calibration in all categories.
Events not assigned to any of the VBF categories are considered for the ggF categories. The two selected large-radius jets are ordered by the value of D bb . The D bb -leading jet is required to have m SD > 50 GeV, while the D bb -subleading jet is required to have m reg > 50 GeV. A BDT is trained to discriminate between the HH signal and QCD multijet or tt background processes. Input variables to the BDT include: the p T and substructure variable τ 32 , used to identify threeprong jets like those arising from fully hadronic top quark decays [86,92], of each of the two selected large-radius jets; the m SD , η, and D bb of the D bb -leading jet; the p miss T ; the p T , η, and mass of the dijet system (m HH ); and the ratios p Tj 1 /m HH , p Tj 2 /m HH , and p Tj 2 /p Tj 1 . The m SD is used as an input to the BDT instead of m reg in order to avoid correlations between the mass regression and the BDT. The simulation modeling of the BDT input variables is validated in signal and background control regions and shown to be accurate. Based on the BDT output score, tight, medium, and loose WPs are defined, corresponding to HH signal selection efficiencies of approximately 23, 27 and 33%, and QCD and tt background misidentification rates of approximately 1, 2 and 12%, respectively. Three SR event categories targeting ggF production are constructed. Events in the highest signal purity category 1 are required to pass the tight BDT WP and the D bb -subleading large-radius jet is required pass the tight D bb WP. Events not in category 1 are grouped into category 2 if they pass the tight BDT WP and the D bb -subleading jet passes the medium D bb WP, or if they fail the tight but pass the medium BDT WP and the D bb -subleading jet passes the tight D bb WP. Finally, events not in categories 1 or 2 are grouped into category 3 if they pass the loose BDT WP and the D bb -subleading jet passes the medium D bb WP. The categorization thresholds are chosen to maximize the expected sensitivity to the SM ggF HH signal.
After all selection requirements, the remaining SM background mainly consists of two processes: tt and QCD multijet production. These background contributions are estimated by fitting the corresponding distributions simultaneously via a maximum likelihood fit to data as detailed in the following, similarly to previous CMS searches for boosted Higgs boson signals [93-95]. Other backgrounds include single Higgs boson production in the ttH and VH modes, V+jets production, and diboson production. They are estimated from simulation in the ggF event categories, and found to be negligible in the VBF event categories.
A control region (CR) enriched in QCD multijet events is selected by changing the requirement on the D bb discriminant. For the VBF categories, this CR is defined by both Higgs boson candidate jets failing the loose D bb WP. Since the D bb discriminant is designed and measured to be independent of the jet mass, the sidebands of the mass distribution of the p T -subleading large-radius jet, namely m reg ∈ [50, 100] ∪ [145, 200] GeV, are used to obtain transfer factors that account for the ratio of selection efficiencies in the QCD-enriched low-D bb CR and the high-D bb VBF SRs. A separate transfer factor is derived using the p T -subleading large-radius jet mass sidebands for each m HH interval in each search category (HP, MP, and LP). For the ggF categories, the QCD-enriched CR is defined by the D bb -subleading large-radius jet failing the medium D bb WP and passing the loose BDT WP. The QCD multijet background in the ggF SRs is estimated using the m reg distribution of the D bb -subleading jet, assuming a constant transfer factor [95]. We have checked whether a higher-order polynomial dependence on the mass shape is necessary using a Fisher F-test [96], and Kolmogorov-Smirnov [97] and saturated model [98,99] goodness-of-fit tests, but a constant factor is found to be sufficient.
For the VBF categories, an auxiliary sample enriched in tt events containing one leptonically decaying W boson is used to extract two corrections for the tt background estimate: one for the difference in H → bb misidentification efficiency in tt events between data and simulation, measured separately for each D bb WP, and the other for the overall normalization of the tt process, measured inclusively. We define this sample with a set of selections closely following the definition of the tt-enriched region in Ref. [86]. For the ggF categories, the simulation predictions of the BDT and m reg distributions for the top quark background are validated and corrected using a similar tt sample containing one leptonic W boson decay, with τ 32 < 0.46 required for two large-radius jets.
The dominant uncertainty in the analysis, corresponding to 62% relative to the extracted signal yield is the statistical uncertainty due to limited size of the data sample. Other sources of uncertainty include: the QCD multijet background estimate (30%); the simulation modeling of the D bb shape and selection efficiency (15%); the jet energy and mass scale and resolution (13%); the theoretical uncertainties from the PDFs, missing higher-order QCD corrections, and initial-and final-state radiation (11%); the modeling of the top quark background (7%); the trigger efficiency measurement (4%); the integrated luminosity measurement [100][101][102] (2%); and the pileup modeling (2%). The impact of each uncertainty is evaluated by computing the uncertainty excluding that source and subtracting it in quadrature from the total uncertainty.
A binned maximum likelihood fit to the observed m HH , D bb , m reg , and BDT distributions is performed using the sum of the signal and background contributions. Because part of the VBF signal can enter the ggF categories, e.g., caused by failing to reconstruct the VBF jets or the Higgs boson candidate jet p T being too low, the fit is performed simultaneously in all ggF and VBF categories, including the CRs. The results of the fit in the ggF BDT event category 1, which dominates the overall sensitivity, are shown in Fig. 1. The results of the fit in the VBF LP, MP, and HP categories are shown in Fig. 2. 138 fb Figure 1: The data and fitted signal and background distributions for the D bb -subleading jet regressed mass are shown for the ggF BDT event category 1, the category accounting for most of the sensitivity to the ggF HH signal. The SM HH (κ 2V = κ V = κ λ = 1) signal is shown scaled to the best fit signal strength µ = 3.5. The lower panel shows the ratio of the data and the total prediction, with its uncertainty represented by the shaded band. The error bars on the data points represent the statistical uncertainties.
The test statistic chosen to determine the signal yield is based on the profile likelihood ratio [103]. Systematic uncertainties are incorporated into the analysis via nuisance parameters and treated according to the frequentist paradigm. To ensure that the background estimation methods and the resulting statistical model are able to disentangle possible signal from background, we perform signal injection tests. We generate various pseudo-experiment data sets corresponding to different signal hypotheses, fit these data sets, and verify that the extracted signal parameters correspond to the injected ones within the fitted uncertainties. Then we pro-  ceed to fit the observed data, for which the signal strength modifier, defined as the ratio of the measured to the expected SM HH production, extracted from a combined fit of all categories is found to be µ = 3.5 +3.3 −2.5 . The best fit value of each parameter of interest and approximate 1, 2, 3, and 5 standard deviation (σ) confidence level (CL) intervals are extracted following the procedure described in Section 3.2 of Ref. [104]. Figure 3 shows the profile likelihood test statistic scan in data over the κ λ -κ 2V plane.
Upper limits on the HH production cross section at 95% CL based on the CL s criterion [105,106] are obtained using asymptotic formulas [107]. The upper limit on the total HH production cross section is observed (expected) to be 9.9 (5.1) relative to the SM prediction. The observed (expected) 95% CL interval for κ λ is found to be [−9.9, 16 [108].
In summary, a search for nonresonant Higgs boson pair (HH) production via gluon-gluon fusion and vector boson (V) fusion in the final state with two bottom quark-antiquark (bb) pairs has been presented. The search is focused on the phase space region where both Higgs bosons are highly Lorentz boosted so that each Higgs boson is reconstructed as a large-radius jet. A novel algorithm based on a graph neural network is applied to identify the jets that correspond to H → bb decays. The data are found to agree with the background-only hypothesis, and an observed (expected) upper limit at 95% confidence level is set to 9.9 (5.1) relative to the standard model cross section. This represents a factor of 30 improvement over the previous best search for a pair of boosted H → bb jets, which used only data collected in 2016 (36 fb −1 ) and less advanced methods for H → bb identification and event selection [94]. Upper limits on the production cross section are set as a function of the coupling modifier parameters κ λ and κ 2V , which parametrize the strengths of the Higgs boson self-coupling, and the quartic VVHH couplings, respectively, relative to their standard model values. The values of κ λ and κ 2V are observed (expected) to be in the ranges [−9.9, 16 Figure 3: Two-parameter profile likelihood test statistic (−2∆ ln L) scan in data as a function of κ λ and κ 2V . The black cross indicates the minimum, while the red diamond marks the SM expectation (κ λ = κ 2V = 1). The gray solid, dashed, dotted, and dash-dotted contours enclose the 1, 2, 3, and 5 σ CL regions, respectively.

fb
Grid and other centers for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies:        [38] CMS Collaboration, "Performance of the CMS muon detector and muon reconstruction with proton-proton collisions at √ s = 13 TeV", JINST 13 (2018) P06015, doi:10.1088/1748-0221/13/06/P06015, arXiv:1804.04528.