Combination of searches for Higgs boson pair production in proton-proton collisions at $\sqrt{s} =$ 13 TeV

This Letter describes a search for Higgs boson pair production using the combined results from four final states: bb$\gamma\gamma$, bb$\tau\tau$, bbbb, and bbVV, where V represents a W or Z boson. The search is performed using data collected in 2016 by the CMS experiment from LHC proton-proton collisions at $\sqrt{s} =$ 13 TeV, corresponding to an integrated luminosity of 35.9 fb$^{-1}$. Limits are set on the Higgs boson pair production cross section. A 95% confidence level observed (expected) upper limit on the nonresonant production cross section is set at 22.2 (12.8) times the standard model value. A search for narrow resonances decaying to Higgs boson pairs is also performed in the mass range 250-3000 GeV. No evidence for a signal is observed, and upper limits are set on the resonance production cross section.


1
The discovery of the Higgs boson (H) by the ATLAS and CMS Collaborations [1-3] was a major step in the understanding of the mechanism of electroweak symmetry breaking. The Higgs boson mass was jointly determined by the two experiments using the CERN LHC Run 1 data to be m H = 125.09 ± 0.24 GeV [4] and recently by CMS using partial Run 2 data with even better accuracy, m H = 125.26 ± 0.21 GeV [5]. With the Higgs boson mass known with a precision better than 0.2%, the structure of the Higgs scalar field potential and the Higgs boson selfcouplings are precisely predicted in the standard model (SM). The current measurements of the properties of the Higgs boson are compatible with the SM predictions [6,7]. The measurement of the Higgs boson self-coupling provides an independent test of the SM and verification that the Higgs mechanism is truly responsible for the electroweak symmetry breaking by giving access to the shape of the Higgs scalar field potential [8]. Access to the Higgs boson trilinear coupling can be obtained by measuring the production of pairs of Higgs bosons (HH) at the LHC. At the same time, several theories predict heavy resonances that decay into HH [9][10][11][12][13][14][15][16]. Studies of pair production of Higgs bosons, each of which can decay in different channels, allow one to probe different regions of the anomalous couplings space and of the resonant invariant mass spectrum. A combination of different channels is therefore needed to obtain the best possible sensitivity for the HH production.
Physics beyond the SM (BSM) can significantly modify the cross section and the kinematic properties of nonresonant Higgs boson pair production. In order to provide model independent constraints on these effects, we introduce an effective field theory (EFT) Lagrangian that extends the SM with dimension-6 operators [25]. This approach results in five anomalous Higgs boson couplings relevant for HH production: the H coupling to the top quark, y t ; the trilinear coupling, λ HHH ; and three additional couplings, denoted by c 2 ; c 2g ; and c g , that represent, respectively, the interactions of a top quark pair with a Higgs boson pair, of a gluon pair with a Higgs boson pair, and of a gluon pair with a single H [17]. We define k λ = λ HHH /λ SM and k t = y t /y SM t . Since a full five-dimensional scan of all possible coupling combinations would be computationally excessive, a clustering strategy [26] has been developed to group together possible combinations of coupling values that present similar kinematic properties. Twelve clusters have been identified, in addition to the SM and the λ HHH = 0 scenarios. Within each cluster, a representative point in the EFT space that we refer to as a benchmark is selected. Each benchmark thus represents a possible modification of the HH signal yield and kinematic distributions due to BSM effects.
In the searches discussed in the Letter, the resonant signal is represented by either a CP-even particle of spin-0 (radion) or spin-2 (graviton) with a width that is much smaller than the detector resolution for the whole range under study.
The ATLAS and CMS Collaborations performed studies of Higgs boson pair production at √ s = 8 TeV [27][28][29]. Limits on nonresonant HH production were set by the CMS Collaboration in the bbγγ and bbττ final states. Here and in the rest of the text the indication of the charge of the decay products is omitted for simplicity of notation. The combination of those searches allowed the CMS Collaboration to set an upper limit on HH production at 43  The CMS apparatus features 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, and a brass and scintillator hadron calorimeter, 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 gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. Events of interest are selected using a two-tiered trigger system [40]. The particle-flow algorithm [41] aims to reconstruct and identify each individual particle in an event with an optimized combination of information from the various elements of the CMS detector. Dedicated b tagging algorithms [42] are used to identify jets originating from b quarks (b jets). 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. [43].
With the current data, the bbγγ [30] analysis is the most sensitive to SM Higgs boson pair production in CMS. Despite a low branching fraction (B(HH → bbγγ) = 0.26%), the analysis profits from a very small background contribution, thanks to the excellent diphoton mass resolution of the CMS experiment for this channel (≈1.6 GeV [30]). To exploit this feature, the analysis relies on a 2D fit of the H → γγ and H → bb invariant mass distributions, where the background coming from the Nγ + jets (N = {0, 1, 2, ..}) continuum is estimated from the mass sidebands. Other background contributions can arise from single H production and constitute up to a third of the total background in the nonresonant searches. They are modeled from Monte Carlo (MC) simulations, including gluon fusion, vector boson fusion, VH, ttH, and bbH production processes. In the nonresonant searches events are classified into low-and highmass categories according to the HH pair reduced massM X = M jjγγ − M jj − M γγ + 250 GeV, which increases the sensitivity to HH production. A further categorization based on the purity of the events is applied to both resonant and nonresonant searches. The signal purity is estimated by means of a boosted decision tree (BDT) built from the b tagging probability of each jet, the angles between the products of the HH system decay, and the transverse momentum of each H candidate. The category most sensitive to nonresonant production is the one with the highest signal purity and reduced mass, with a SM signal over background ratio (S/B) of ≈5%. The low reduced mass categories enhance sensitivity to large k λ values. In the resonant search, two different BDTs are trained for resonance masses m X higher or lower than 600 GeV.
The bbττ analysis [32] combines a relatively low background contamination with a relatively large branching fraction (7.3%), for a final S/B of ≈0.4% in the most sensitive category. At least one isolated hadronically decaying τ lepton must be present in the event, together with a second isolated lepton that is oppositely charged. The second lepton can be either an electron or a muon (semileptonic final state) or another hadronically decaying τ lepton (fully hadronic final state). Events are categorized according to the number of b-tagged jets (one or two) in the event. Events with a boosted H jet candidate are assigned to a dedicated boosted category. A BDT discriminant based on the kinematic differences between the HH and tt processes is used in the semileptonic channels in order to reduce the large amount of tt background. A kinematic fit [44] is performed to reconstruct the Higgs boson mass. For the nonresonant searches, the stransverse mass m T2 [32, 45,46] is found to provide the best separation between the background and the signal [45]. For the resonant searches, a further kinematic fit, detailed in Ref. [47], is performed to reconstruct the most probable mass of the resonance. The largest background contributions come from Drell-Yan, tt, processes and events comprised uniquely of jets produced via the strong interaction (QCD multijet). Additional background sources considered are VH, VV, and W+jets. The Drell-Yan background is estimated from a leading order simulation with a scale factor obtained from data in a Z boson enriched control region to account for higher-order corrections. The QCD multijet background is estimated using three control regions with different lepton charge (same-sign region) and isolation requirements. All other background contributions are estimated from MC.
The bbbb final state has the highest HH decay branching fraction (33.6%). Three optimized resonance searches are performed targeting resonances of different masses. For a resonance mass below 0.7 TeV, four b-tagged jets [34] are required in the final state. For high-mass searches (m X > 1.2 TeV), selected events must contain two jets, each associated with a boosted H decaying into a bb pair merged in a single jet ("H jet candidate") [35]. This candidate reconstruction uses jet substructure and jet flavor tagging techniques [48]. The search sensitivity in the intermediate mass range, 0.7-1.2 TeV, is improved by considering both the final state with four b jets and the case with one Higgs jet candidate and two b jets [36]. A dedicated nonresonant HH search is performed in the four b jet final state [37]. Triggers using b tagging and jet substructure techniques are used to collect these events. The sensitivity of the analysis is improved using a BDT that employs jet-related and HH decay kinematic variables and global event properties. Under SM assumptions, the S/B can be as high as 1% in the most sensitive BDT bins. Further sensitivity to BSM nonresonant HH production, which often results in boosted topologies, is obtained adding the final states with one or two H jet candidates [36]. The dominant background comes from QCD multijet events and it is estimated using data sideband regions [34-36] and a hemisphere mixing technique [37]. The residual background is dominated by tt events and is estimated using simulation. Such events correspond to 10, 15, and 1% of the total background in the low-, intermediate-, and high-mass ranges, respectively, for resonant searches and 5 (15)% for nonresonant searches in the final states with four b-tagged jets (with H jet candidates).
The bbVV analysis [39] includes the HH → bbWW → bb ν ν and the HH → bbZZ → bb νν processes, for a total branching fraction B(HH → bb ν ν) = 2.72% when all possible lepton flavor decays are considered (including the τ leptonic decays). This channel has large backgrounds, which are predominantly from tt and Drell-Yan processes. All background sources are estimated from MC with the exception of Drell-Yan, which is estimated from data. In order to improve the rejection of these large backgrounds, a deep neural network (DNN) discriminant technique has been implemented. Two separate parametrized DNNs have been trained, one for the resonant and one for the nonresonant search. The first is parametrized according to the mass of the resonance, while the second depends on the parameters k λ and k t . For each nonresonant benchmark, the closest point in the (k λ , k t ) plane, according to the statistical measure defined in Ref. [26], has been used to define the DNN parameters.
For both the resonant and nonresonant searches, likelihood fits are performed using the statistical toolkits, ROOFIT [49] and ROOSTATS [50]. The signal strength (µ), defined as the ratio between the observed and expected signal rates, is estimated with its corresponding confidence interval via the profile likelihood ratio test statistic [51]. The latter depends on the signal strength as well as on the nuisance parameters, which account for various experimental and theoretical uncertainties. The reference value for the expected signal strength is chosen as the SM gluon fusion HH production cross section (33.53 fb) in nonresonant searches. The likelihood fits are performed with respect to the observed data or the Asimov data set [52] for assessing expected results using the asymptotic approximation [52,53]. Limits are set at 95% confidence level (CL) using the CL s criterion [54,55]. For all measurements, the H mass is fixed at m H = 125 GeV, and branching fractions are assumed to be equal to the SM predictions. When investigating signal models corresponding to the shape benchmarks, the single H production cross sections are all assumed to have their SM values.
When combining results, the systematic uncertainties from various sources are accounted for as follows. A polynomial interpolation between alternative shape variations is used to model systematic effects on the shape of the discriminant variables. Lepton and photon reconstruction and identification efficiencies and energy scale corrections are assumed to be fully correlated across the channels and analysis categories that use the same objects. The uncertainties in the jet energy scale corrections are divided into multiple sources. The effect of each source on the rate and shape of the final observable for different processes is considered in the bbVV analysis and is treated as fully correlated with the normalization effects in the bbττ and bbγγ channels. In these latter channels, the shape effects from the individual sources have been verified to be negligible, and only the cumulative effects on the shape are considered. They are assumed to be correlated with the normalization and shape variation from the bbbb analysis. Jet energy resolution effects are negligible in the bbττ channel and assumed to be fully correlated between the bbVV and the bbbb channels. The effects of jet energy scale and resolution are included in the functional form for the shape used in the bbγγ channel and are assumed to be uncorrelated with the other channels. Uncertainties related to the b tagging are dominated by the modeling of heavy flavor production when measuring data/MC scale factors, and are considered correlated across the bbττ, bbγγ, resonant bbbb, and bbVV channels. The nonresonant bbbb analysis uses a different b tagging method and is considered uncorrelated.
An integrated luminosity uncertainty of 2.5% [56] is applied in a fully correlated way to all channels and all processes estimated from simulation. Different uncertainties are applied to background processes that are estimated from data. The uncertainties in the total cross sections of the common background processes are assumed to be correlated across all channels, and are of the order of 5% for the most relevant ones (single top, tt, VH). The uncertainty in the same-sign to opposite-sign candidate ratio in the bbττ channel is propagated to the estimation of the multijet background, and the uncertainties in the scale factors applied to the Z/γ * → background estimation to correct for higher-order effects are also taken into account. A normalization uncertainty is considered in the bbbb nonresonant search to account for residual biases in the hemisphere mixing technique. Some background estimates in the bbττ, bbVV, and nonresonant bbbb analyses have nonnegligible statistical uncertainties due to the limited number of events passing the selection in data control regions or simulated data samples. These are taken into account by allowing each bin in each template shape to fluctuate independently according to a Poisson distribution. These uncertainties are assumed to be uncorrelated across bins in the individual template shapes. The nonresonant signal uncertainties include contributions coming from variations in the renormalization and factorization scales, these amount to +4.3% −6.0% [17,20] of the nonresonant signal cross section. Other theoretical uncertainties such as those in α S , PDFs, and finite top quark mass effects at next-to-next-to-leading order result in a further 5.9% uncertainty [18,19,21]. These effects are assumed to be fully correlated across the different channels. The α S , PDFs, and scale variation effects are also included for single H background contributions. These uncertainties are considered fully correlated for the VH production in the bbττ, bbVV, and bbγγ channels. The uncertainty in the branching fraction of the Higgs boson to bb [17] is also assumed to be fully correlated across all channels.
The event yield in data is small, so the statistical uncertainties are much larger than the systematic ones. Those   the yield in the most sensitive bins of the BDT and the overall background normalization in the bbbb channel; the hadronically decaying τ lepton energy scale effects in the bbττ analysis; and the uncertainty in the functional form used to model the signal shape in the bbγγ channel. These effects are as large as 10 (5)% for the bbbb and bbττ (bbγγ) uncertainties. Due to its lower overall sensitivity, the systematic uncertainties affecting the bbVV analysis have little effect on the combined result. The largest sources of systematic uncertainty for this channel arise from the uncertainties in the tt cross section, electron identification efficiency, and b tagging efficiency.
With all the correlations across the channels included, the observed and expected limits at 95% confidence level on the nonresonant HH production signal strength are measured to be 22.2 and 12.8 times the SM expectations, respectively. They are shown in Fig. 1 for the individual channels and their combination. Small excesses, compatible with statistical fluctuations, are observed in the bbbb, bbττ, and bbγγ final states and result in a small excess in the combined result. A scan is performed for different values of the k λ parameter, while keeping all other EFT parameters fixed at their SM values. The value of k λ affects both the expected cross section and the HH decay kinematics. For each value, these differences are fully simulated and considered in the scan. Resulting limits are reported in Fig. 2. The exclusion limit as a function of k λ closely follows the features of the HH production cross section and HH invariant mass distribution M H H [57], which are sculpted by the interference between the HH production via the trilinear Higgs coupling and the emission of an HH pair from a top quark loop. The minimum at k λ = 2.46 corresponds to the maximum negative interference between the two diagrams, which results in a minimum of the cross section but at the same time enhances the relative importance of tails in the M H H distribution. The maximum at k λ ≈5 is due to the softness of the M H H spectrum for such values of the trilinear coupling, causing a larger fraction of events to fall outside experimental acceptance. As |k λ | increases, the production via the trilinear Higgs coupling becomes dominant and the limit asymptotically approaches the same value for both k λ −10 and k λ 10. This is reflected in the observed exclusion limit as well, where the significance of the small observed excess is relatively less important in the more sensitive small k λ region than at large values of k λ . When fixing all the other EFT parameters to their SM values, the k λ parameter is observed (expected) to be constrained to the range −11.8 < k λ < 18.8 (−7.1 < k λ < 13.6) at 95% CL. The observed exclusions for the different EFT benchmarks [26] are in the range of 100-3000 fb, and can be seen in Appendix A. A small excess, similar to that observed at the SM value, is present across most of the phase space with the exception of the more boosted topologies.
The resonant search is performed in the range of masses from 250 to 3000 GeV. Under the hypothesis of a narrow-width resonance, no significant excess is found across the whole range for either a spin-0 or a spin-2 resonance. The results of the combined resonant search are shown in Fig. 3 for the spin-0 model, and in Appendix A for the spin-2 case.
In summary, a combination of searches for nonresonant and resonant Higgs boson pair production has been presented. The combination includes the bbγγ, bbττ, bbbb, and bbVV channels, where V represents a W or Z boson, using a data sample collected in proton-proton collisions at √ s = 13 TeV, which corresponds to an integrated luminosity of 35.9 fb −1 . Upper limits at 95% confidence level (CL) on the Higgs boson pair production cross section are obtained. For the nonresonant production mechanism, the observed (expected) 95% CL corresponds to 22.2 (12.8) times the theoretical prediction for the standard model cross section. An effective field theory framework is used to parametrize the cross section as a function of anomalous couplings of the Higgs boson. When varying only the ratio k λ between the Higgs trilinear coupling λ HHH and its standard model (SM) expectation, values of k λ in the region −11.8 < k λ < 18.8 (−7.1 < k λ < 13.6) are still allowed by the observed (expected) data. For the resonant production mechanism, upper exclusion limits at 95% CL are obtained for the production of a narrow resonance with mass ranging from 250 to 3000 GeV. These results represent both the most sensitive and most comprehensive study of double Higgs production at the LHC to date.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. [6] ATLAS and CMS Collaborations, "Measurements of the Higgs boson production and decay rates and constraints on its couplings from a combined ATLAS and CMS analysis of the LHC pp collision data at √ s = 7 and 8 TeV", JHEP 08 (2016) 045, doi:10.1007/JHEP08(2016)045, arXiv:1606.02266.
[31] ATLAS Collaboration, "Search for Higgs boson pair production in the γγbb final state with 13 TeV pp collision data collected by the ATLAS experiment", (2018). arXiv:1807.04873.
[33] ATLAS Collaboration, "A search for resonant and non-resonant Higgs boson pair production in the bbτ + τ − decay channel in pp collisions at √ s = 13 TeV with the ATLAS detector", (2018). arXiv:1808.00336. Submitted to Phys. Rev. Lett. [36] CMS Collaboration, "Search for production of Higgs boson pairs in the four b quark final state using large-area jets in proton-proton collisions at √ s = 13 TeV", (2018). arXiv:1808.01473. Submitted to JHEP.
[38] ATLAS Collaboration, "Search for pair production of Higgs bosons in the bbbb final state using proton-proton collisions at √ s = 13 TeV with the ATLAS detector", (2018). arXiv:1804.06174. [40] CMS Collaboration, "The CMS trigger system", JINST 12 (2017) P01020, doi:10.1088/1748-0221/12/01/P01020, arXiv:1609.02366. Figure 4 shows the 95% CL upper limit on the HH production cross section for each benchmark topology, the standard model and the k λ = 1 case and for each final state. As can be seen in the plot, different benchmark shapes result in very different sensitivities for a given analysis. The combined limits alone are shown in Fig. 5.

A Supplemental material
The sensitivity of each analysis to the SM HH production is shown in Fig. 6, together with the 8 TeV CMS result.
The exclusion limits for a spin-2 narrow width resonance are shown in Fig. 7.  Figure 4: The 95% CL exclusion limits on the nonresonant Higgs boson pair production cross sections for different EFT benchmark topolo-gies (bins 1 to 12) as defined in Ref. [26]. Each benchmark represents a possible modification in the HH signal yield and kinematic distribu-tions due to BSM effects. The last two bins show the 95% CL exclusion limits for the SM and for the k λ = 0 scenarios. Limits are shown for each of the four final states separately and for the combination. Each bin of the histogram shows the limit for a different benchmark (bins 1 to 12) as defined in Ref. [26], for SM production and for k λ = 0 production (last 2 bins). Each benchmark represents a possible modification in the HH signal yield and kinematic distributions due to BSM effects.