Combination of searches for Higgs boson pair production in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector

This Letter presents results from a combination of searches for Higgs boson pair production using 126$-$140 fb$^{-1}$ of proton-proton collision data at $\sqrt{s}=13$ TeV recorded with the ATLAS detector. At 95% confidence level (CL), the upper limit on the production rate is 2.9 times the standard model (SM) prediction, with an expected limit of 2.4 assuming no Higgs boson pair production. Constraints on the Higgs boson self-coupling modifier $\kappa_{\lambda}=\lambda_{HHH}/\lambda_{HHH}^\mathrm{SM}$, and the quartic $HHVV$ coupling modifier $\kappa_{2V}=g_{HHVV}/g_{HHVV}^\mathrm{SM}$, are derived individually, fixing the other parameter to its SM value. The observed 95% CL intervals are $-1.2<\kappa_{\lambda}<7.2$ and $0.6<\kappa_{2V}<1.5$, respectively, while the expected intervals are $-1.6<\kappa_{\lambda}<7.2$ and $0.4<\kappa_{2V}<1.6$ in the SM case. Constraints obtained for several interaction parameters within Higgs effective field theory are the strongest to date, offering insights into potential deviations from SM predictions.

Since the discovery of the Higgs boson () by the ATLAS and CMS Collaborations at the Large Hadron Collider (LHC) in 2012 [1,2], understanding its intrinsic properties and interactions has been a priority.The Higgs self-coupling is directly related to the shape of the Higgs scalar field potential, which is important for understanding the mechanism of electroweak symmetry breaking and serves as an essential test of the electroweak theory.After the symmetry breaking, the Higgs potential predicted in the standard model (SM) can be expanded in the Higgs boson field  near its minimum:  () = 1 2  2   2 +      3 + O ( 4 ), where   is the Higgs boson mass and  ≈ 246 GeV is the field's vacuum expectation value [3].The Higgs boson's trilinear self-coupling  SM    is equal to  2  /2 2 , its coupling  SM  to vector bosons ( = , ) is equal to 2 2   /, and its coupling  SM   f to fermions is equal to   /.The quartic coupling between two Higgs bosons and two vector bosons,  SM   , is equal to  SM  / [3].The production of Higgs boson pairs () via gluon-gluon fusion (ggF) and vector-boson fusion (VBF) provides a direct probe of     and    , which affects the pair-production differential cross section at tree level.Observed (expected) 95% confidence level (CL) upper limits on the  production cross section have been set at 2.4 (2.9) and 3.4 (2.5) times the SM prediction by previous ATLAS [4] and CMS [5] search combinations, respectively.
Deviations from the SM can be expressed in terms of coupling modifiers   [6,7] or Wilson coefficients   in the Higgs effective field theory (HEFT) [8,9], as illustrated in Fig. 1.In the  framework, the coupling modifiers   ,   ,   , and  2 are each defined as the ratio of the Higgs boson coupling to its SM value, e.g.,   =     / SM    .In the HEFT framework, new physics in the electroweak sector is described through anomalous couplings of the Higgs boson.The organization of the HEFT Lagrangian is guided by chiral perturbation theory, with the low-energy dynamics of electroweak symmetry breaking described using a nonlinear realization of the gauge symmetry group SU(2)  × U(1)  .One advantage of the HEFT framework is that the anomalous single-Higgs-boson and  couplings are defined separately, allowing simplified  interpretations.In the HEFT Lagrangian, ggF  production is described at leading order by five relevant operators, and their associated Wilson coefficients are   ℎ ,  ℎ ,  ℎℎℎ ,  ℎℎ , and   ℎℎ .In this formalism,   ℎ and  ℎℎℎ are equivalent to   and   .In this analysis,   ,   (  ℎ ), and  ℎ are set equal to the SM predictions, because those parameters are constrained by precise measurements of single-Higgs-boson production [5,10].
In the SM, destructive interference between the ggF  production diagrams in Figs.1(a) and 1(b) makes softer Higgs bosons more sensitive for constraining   .For   = 125 GeV and √  = 13 TeV proton-proton collisions, the predicted cross section is  SM ggF () = 31.1 +1.9 −7.1 (scale +  top ) ± 0.9 (PDF +  s ) fb [11][12][13][14][15][16][17][18] at next-to-next-to-leading order in  s and including and including an uncertainty related to the choice of the virtual top-quark mass scheme [18].The "PDF +  s " uncertainty accounts for uncertainties in the parton distribution functions and strong coupling constant, the "scale" uncertainty is due to the finite order of quantum chromodynamics (QCD) calculations, and the " top " uncertainty is related to the top-quark mass scheme.For SM VBF  production, divergences in the diagrams shown in Figs.1(d) and 1(e) cancel out due to perturbative unitarity.If  2 deviates from the SM prediction, this cancellation no longer occurs, leading to a linear dependence of the cross section on the effective center-of-mass energy of the incoming vector bosons [19].Consequently, the Higgs bosons are expected to be more energetic in non-SM scenarios.The cross section for VBF  production is  SM VBF () = 1.73 ± 0.04 fb at next-to-next-to-next-to-leading order in QCD [20][21][22][23][24].This Letter presents a combination of results from the  b b [25, 26],  b +  − [27],  b [28], multilepton [29], and  bℓℓ +  miss T [30] decay channels, probing more than half of the  decays.The first three analyses have been improved since the previous combination [4], and the other two are newly included.The  →  b b decay mode has the advantage of having the largest SM  decay branching fraction (33.9%), but it also has the largest SM background, due to the abundance of QCD multĳet events.Given its capability to probe relatively high-energy Higgs bosons, both the resolved [25] and boosted topologies [26] are now used to reconstruct the Higgs bosons.The  →  b +  − decay mode has one of the larger branching fractions (7.3%) among the investigated  decay channels and benefits from having only moderate background contamination.In the corresponding search [27], one of the  leptons is required to decay hadronically, ensuring orthogonality with the  bℓℓ +  miss T search.Although the  →  b decay mode has a small branching fraction (0.26%), it has high trigger efficiency and a clean experimental signature.The  b +  − [27] and  b [28] analyses have been improved through optimized classification of selected events to enhance the sensitivity to the Higgs boson couplings.Furthermore, the  b +  − analysis now benefits from more accurate background modeling and larger samples of simulated events.The multilepton analysis is designed to select  events in  b  * ,  *  * ( =  or ),  *  +  − ,  +  −  +  − ,  * , and  +  − decay channels with leptons in the final states; the total branching fraction is around 6.5%.The  bℓℓ +  miss T search targets final states arising from  decay channels where one of the Higgs bosons decays to a -quark pair and the other to either a boson pair (  * ,  * ) or a -lepton pair, which then decays to a pair of opposite-sign leptons (ℓ = , ) and neutrinos, for a total branching fraction of 2.9%.Depending on the analysis, the final discriminating variable can be the  invariant mass, the diphoton invariant mass, or the multivariate classifiers used to separate signal from background.
The analyses under consideration use the full sample of √  = 13 TeV proton-proton ( ) collision data recorded with the ATLAS detector during run 2 of the LHC.The integrated luminosity ranges from 126 to 140 fb −1 depending on the trigger selection [31].The ATLAS experiment is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and nearly 4 coverage in solid angle [32][33][34].A software suite [35] is used in data simulation, in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.The searches use a common set of event generators to describe ggF and VBF  production in the   collisions.Reweighting methods are used to estimate the total and differential signal yields at a given value of   from samples simulated for different values of   and  2 [4] or to estimate the particle-level    distributions for alternative values of the Wilson coefficients using parameters from Ref. [36].
The results are derived from a likelihood function (, ), where  denotes the vector of parameters of interest (POIs) in the statistical model and  is a set of nuisance parameters (NPs), including systematic uncertainty contributions and background parameters.This global likelihood function is the product of individual search likelihoods.The profile-likelihood-ratio test statistic is used to determine the 68% and 95% CL intervals and local significance in the asymptotic approximation [37].
The CL s method [38] is utilized to derive upper limits on the  production cross section.To evaluate the expected limits, Asimov datasets [37] are generated, setting all NPs to their best-fit values in data and fixing the POIs to those posited in the hypothesis under test.The event samples from the combined searches are scrutinized for overlaps in both real and simulated data; they are found to be less than 1% in the signal regions, and, thus, considered negligible.
Complete discussions of the systematic uncertainties considered in the individual searches are provided in Refs.[25][26][27][28][29][30].Correlations of these uncertainties between different searches are investigated.Uncertainties related to the data-taking conditions, such as those associated with the integrated luminosity and the mismodeling of the multiple   interactions per bunch crossing, are assumed to be correlated across the searches.An exception is the integrated luminosity uncertainty in the resolved  b b analysis [25], which employs a different calibration version.Where applicable, uncertainties associated with physics objects common to two or more searches are considered correlated.Correlations are also assumed for theoretical uncertainties affecting simulated signal and background processes, such as uncertainties in the QCD scale, proton parton distribution functions, and Higgs boson decay branching fractions.Systematic uncertainties that significantly influence the individual searches but are strongly constrained or pulled in the data fitting are treated as uncorrelated to prevent undue influence on the other searches.However, the impact of treating them as correlated or uncorrelated in the combination was checked and found to be negligible.
The signal strength    is defined as the ratio of the measured inclusive ggF and VBF  production cross section to the SM prediction  SM ggF + VBF () = 32.8+2.1 −7.2 fb.This    measure assumes that the relative ggF and VBF production cross sections, Higgs boson decay branching fractions, and relative kinematic distributions correspond to the SM predictions.The fit to data indicates a value of    = 0.5 +1.2 −1.0 = 0.5 +0.9 −0.8 (stat) +0.7 −0.6 (syst), where "stat" and "syst" denote the statistical and systematic uncertainties, respectively.The result is compatible with the SM prediction, with a  value of 0.64.Assuming  ggF + VBF () =  SM ggF + VBF (), the expected value is    = 1.0 +1.2 −1.0 = 1.0 +1.0 −0.9 (stat) +0.7 −0.5 (syst).The primary systematic uncertainty arises from an estimated uncertainty of 100% in modeling the radiation of additional heavy-flavor jets in the ggF single-Higgs-boson background production process [39][40][41][42][43], affecting    by 25%.The observed (expected) significance of    is 0.4 (1.0) standard deviations, with respect to the hypothesis of no  production.No significant  signal is observed above the expected background, and a 95% CL upper limit of 2.9 is placed on    .If  production is absent, the expected Obs.Exp.
Observed limit (95% CL) Expected limit (95% CL) ( HH = 0 hypothesis) Expected limit ±1 Expected limit ±2 Figure 2: Observed and expected 95% CL upper limits on the signal strength for inclusive ggF  and VBF  production from the  b +  − ,  b,  b b, multilepton and  bℓℓ +  miss T decay channels and their statistical combination.The predicted SM cross section assumes   = 125 GeV.The expected limit, along with its associated ±1 and ±2 bands, is calculated for the assumption of no  production and with all NPs profiled to the observed data.
95% CL upper limit is 2.4, and in the SM case (   = 1) the expected upper limit is 3.4.The expected upper limit is 17% lower than in the previous combination [4]: 13% from improvements in the  b +  − ,  b and  b b analyses and an additional 4% from the inclusion of the multilepton and  bℓℓ +  miss T channels.This combination provides the best expected sensitivity to the  production cross section to date. Figure 2 displays the limits from the individual searches and their combination [44] highlighting the  b +  − channel as the one expected to constrain    the most.The  value for compatibility between the    value measured in the combination and those measured in the individual searches is 0.16.The observed and expected 95% CL upper limits on  ggF + VBF () from the combination are 86 and 71 fb, respectively, derived in this case excluding theoretical uncertainties in the  production cross section.
The self-coupling modifier   is explored in the ggF and VBF  production processes.The impact of   on the single-Higgs-boson background productions and the Higgs decay widths is neglected.Assuming that other Higgs boson couplings conform to the SM predictions, a fit to data yields   = 3.8 +2.1 −3.6 , which is compatible with the SM prediction, with a  value of 0.53.The expected value of   is 1.0 +4.7 −1.5 when assuming SM  production.The observed (expected) 95% CL interval is −1.2 <   < 7.2 (−1.6 <   < 7.2), representing the best expected sensitivity to the Higgs boson self-coupling to date.The values of the test statistic as a function of   are shown in Fig. 3(a) for both the individual searches and their combination, highlighting the  b channel as the most sensitive.Similarly,  2 is explored in the VBF  production process.Assuming the SM predictions for other Higgs boson couplings, the observed (expected) value is  2 = 1.02 +0.22  −0.23 ( 2 = 1.00 +0.40 −0.36 ).The observed (expected) 95% CL interval is 0.6 <  2 < 1.5 (0.4 <  2 < 1.6).The values of the test statistic as a function of  2 are shown in Fig. 3(b), highlighting the  b b analysis as the most sensitive, mainly due to the boosted channel [26].A deficit of data events in this channel results in stronger constraints on  2 than expected.To reduce model dependence, two-dimensional contours of −2 ln Λ in the  2 -  plane are presented in Fig. 3(c).The  value for compatibility of the combined measurement and the SM prediction is 0.78.
For the HEFT interpretation the three most sensitive  decay channels,  b +  − ,  b, and  b b, are combined.The VBF  process is ignored, since it is sensitive only to  ℎℎℎ and the predictions for this Wilson coefficient are not available for this process.One-dimensional constraints are evaluated separately for the coefficients  ℎℎ and   ℎℎ , with all other coefficients fixed to the SM predictions.At 95% CL, the observed interval on  ℎℎ is −0.38 <  ℎℎ < 0.49; if the SM value  ℎℎ = 0 is assumed, the expected 95% CL interval is −0.36 <  ℎℎ < 0.36.Similarly, the observed (expected) 95% CL interval on   ℎℎ is −0.19 <   ℎℎ < 0.70 (−0.27 <   ℎℎ < 0.66).These represent the most stringent constraints to date on  ℎℎ and   ℎℎ .The results are compatible with the SM predictions, with  values of 0.087 and 0.16, respectively.Figure 4 displays the two-dimensional test-statistic contours in the coefficient spaces of ( ℎℎ ,  ℎℎℎ ), (  ℎℎ ,  ℎℎℎ ), and ( ℎℎ ,   ℎℎ ), with each plot fixing   ℎℎ ,  ℎℎ , or  ℎℎℎ , respectively, to its SM value.Two minima are expected because of the quadratic dependence of the cross section on the coefficients.The  values for compatibility of the  ℎℎ - ℎℎℎ ,   ℎℎ - ℎℎℎ , and  ℎℎ -  ℎℎ measurements with the SM predictions are 0.044, 0.21 and 0.031, respectively.The relatively low  values are primarily due to the  b b analysis [25], where the data-driven background modeling cannot perfectly describe the background distribution in data, making non-SM signals more favorable in the fit.Because of insufficient sensitivity, the combination does not allow simultaneous constraints to be placed on  ℎℎℎ ,  ℎℎ , and   ℎℎ in a more model-independent manner [45].
In summary, this Letter presents a combination of the results of searches for  production in the  b +  − ,  b,  b b, multilepton, and  bℓℓ +  miss T decay channels, utilizing the complete LHC run 2 dataset of 13 TeV proton-proton collisions recorded with the ATLAS detector.This new combination provides the best expected sensitivities to the  production cross section and the Higgs boson self-coupling, superseding the results on the di-Higgs measurements of Ref. [4].The results agree well with the SM predictions.When using Higgs effective field theory to interpret the measurements, unprecedented constraints are placed on the effective  and  t interactions.[45] The results in both the  and HEFT frameworks are affected by the issue concerning the ggF  prediction for BSM scenarios in Powheg reported in the erratum to Ref. [47].The main effect is on the high-   region and decreases the signal acceptance.The changes to resolve this issue described in Ref. [48] are not included, to be consistent with individual channels.It was estimated that the impact on the 95% CL limits would be smaller than 4% for   , negligible for  ℎℎ , and at most 10% for   ℎℎ .

Figure 1 :
Figure 1: Leading-order Feynman diagrams showing the production of Higgs boson pairs via the ggF (a), (b), (f)-(h) and VBF (c)-(e) processes.Each diagram is sensitive to specific coupling factors, denoted by   in the  framework or   in the HEFT.Diagrams (a)-(e) occur in the SM predictions, while diagrams (f)-(h) manifest only when deviations from the SM predictions are present in the coefficients  ℎ ,  ℎℎ , or   ℎℎ .
limit on HH signal strength HH Combined

Figure 3 :
Figure 3: Expected values (dashed lines) of the test statistic (−2 ln Λ) as functions of (a)   and (b)  2 .These results are shown for the decay channels  b (purple),  b +  − (green), multilepton (cyan),  b b (blue), and  bℓℓ +  miss T(brown), as well as their combination (black).The observed values from the combined data are depicted by solid black lines.These results are computed with the assumption that all other Higgs boson couplings follow the SM predictions.(c) The expected 95% CL contours in the  2 -  plane, corresponding to the individual decay channels and their combination, are illustrated using dashed lines.The observed contour from the combined results is depicted by a solid black line.The SM prediction is marked by a star, and the combined best-fit value is indicated by a cross.

[ 44 ]
The test statistic and statistical uncertainties of the signal MC samples are updated in the  bℓℓ +  miss T result compared to Ref. [30].
ℎℎ - ℎℎℎ , (b)   ℎℎ - ℎℎℎ , and (c)   ℎℎ - ℎℎ HEFT parameter spaces, with   ℎℎ ,  ℎℎ , and  ℎℎℎ fixed to their SM values, respectively.The corresponding SM expected contours are shown by the inner and outer shaded regions.The SM prediction is indicated by the star, while the best-fit value is shown by the cross.Search for the non-resonant production of Higgs boson pairs via gluon fusion and vector-boson fusion in the b +  − final state in proton-proton collisions at √  = 13 TeV with the ATLAS detector, Phys.Rev. D 110 (2024) 032012, arXiv: 2404.12660[hep-ex].[28] ATLAS Collaboration, Studies of new Higgs boson interactions through nonresonant  production in the  b final state in   collisions at √  = 13 TeV with the ATLAS detector, JHEP 01 (2024) 066, arXiv: 2310.12301[hep-ex].[29] ATLAS Collaboration, Search for non-resonant Higgs boson pair production in final states with leptons, taus, and photons in   collisions at √  = 13 TeV with the ATLAS detector, (2024), arXiv: 2405.20040[hep-ex].[30] ATLAS Collaboration, Search for non-resonant Higgs boson pair production in the 2 + 2ℓ +  miss T final state in   collisions at √  = 13TeV with the ATLAS detector, JHEP 02 (2024) 037, arXiv: 2310.11286[hep-ex].
[25] ATLAS Collaboration, Search for nonresonant pair production of Higgs bosons in the  b b final state in   collisions at √  = 13 TeV with the ATLAS detector, Phys.Rev. D 108 (