Constraints on the Higgs boson self-coupling from single- and double-Higgs production with the ATLAS detector using $pp$ collisions at $\sqrt{s}=13$ TeV

Constraints on the Higgs boson self-coupling are set by combining double-Higgs boson analyses in the $b\bar{b}b\bar{b}$, $b\bar{b}\tau^+\tau^-$ and $b\bar{b} \gamma \gamma$ decay channels with single-Higgs boson analyses targeting the $\gamma \gamma$, $ZZ^*$, $WW^*$, $\tau^+ \tau^-$ and $b\bar{b}$ decay channels. The data used in these analyses were recorded by the ATLAS detector at the LHC in proton$-$proton collisions at $\sqrt{s}=13$ TeV and correspond to an integrated luminosity of 126$-$139 fb$^{-1}$. The combination of the double-Higgs analyses sets an upper limit of $\mu_{HH}<2.4$ at 95% confidence level on the double-Higgs production cross-section normalised to its Standard Model prediction. Combining the single-Higgs and double-Higgs analyses, with the assumption that new physics affects only the Higgs boson self-coupling ($\lambda_{HHH}$), values outside the interval $-0.4<\kappa_{\lambda}=(\lambda_{HHH}/\lambda_{HHH}^{\textrm{SM}})<6.3$ are excluded at 95% confidence level. The combined single-Higgs and double-Higgs analyses provide results with fewer assumptions, by adding in the fit more coupling modifiers introduced to account for the Higgs boson interactions with the other Standard Model particles. In this relaxed scenario, the constraint becomes $-1.4<\kappa_{\lambda}<6.1$ at 95% CL.


Introduction
Since the discovery of the Higgs boson by the ATLAS and CMS collaborations [1,2] at the Large Hadron Collider (LHC) [3], a major goal of the physics programme of the LHC experiments has been to measure its properties and determine whether they correspond to those predicted by the Standard Model (SM) of particle physics [4][5][6][7] or involve new phenomena beyond those described by this theory. One of the most intriguing and interesting characteristics of the SM is that the gauge electroweak (EW) symmetry is broken spontaneously by the non-trivial structure of the Higgs boson [8][9][10][11][12][13] potential, related to its self-interaction. In the SM, this mechanism allows elementary particles to acquire their mass, while preserving perturbative unitarity up to very high energies. The Higgs boson potential also plays a fundamental role in understanding the stability of our universe [14].
The Higgs boson self-interactions are characterised by the trilinear self-coupling . In the SM, the Higgs boson self-coupling can be predicted at lowest order from the values of the Higgs boson mass [15] and the Fermi constant F [16]: At the LHC the Higgs boson self-interaction is directly accessible via the production of Higgs boson pairs (here referred to as double-Higgs production). In this Letter the three most sensitive double-Higgs decay channels,¯,¯+ − , and¯¯ [17][18][19], are combined using the complete dataset collected by ATLAS at √ = 13 TeV in the data-taking period 2015-2018, corresponding to an integrated luminosity of 126-139 fb −1 . This combination is used to place constraints on the double-Higgs production cross-section and on the Higgs boson self-coupling. Results are reported in terms of the coupling modifier defined as the ratio of the Higgs boson self-coupling to its SM value, = / SM .
The Higgs boson self-interaction also contributes to other processes via sizeable next-to-leading-order (NLO) EW corrections. In particular, it has been shown [20-25] that the single Higgs boson (here referred to as single-Higgs) production cross-sections and branching ratios are also modified if the Higgs boson self-coupling deviates from the SM prediction.
More stringent constraints on are also reported in this Letter from combinations of the recent ATLAS single-Higgs results [26] based on the full Run 2 data set from the , * , * , + − and¯decay channels with the above mentioned double-Higgs results. The single-Higgs measurements of the simplified template cross-sections (STXS) and the double-Higgs results have been parameterised to take into account the impact of and the other coupling modifiers. This more comprehensive combination makes it possible to perform tests of relaxing the assumptions about Higgs boson interactions with the other SM particles.
A previous ATLAS combination of searches for non-resonant and resonant pair production was performed on a partial Run 2 dataset, using up to 36.1 fb −1 of data [27]. The combined observed (expected) upper limit on non-resonant production at 95% confidence level (CL) was 6.9 (10) times the predicted SM cross-section. When varying the Higgs boson trilinear self-coupling from its SM value, the allowed range of the self-coupling modifier was observed (expected) to be −5.0 ≤ ≤ 12.0 (−5.8 ≤ ≤ 12.0). The CMS Collaboration also published a combination of searches using its full Run 2 dataset, up to 138 fb −1 of data [28]. The CMS combined observed (expected) upper limit on non-resonant production at 95% CL is 3.4 (2.5) times the predicted Standard Model cross-section, and the observed allowed range of the self-coupling modifier is −1.24 ≤ ≤ 6.49.

Theoretical framework
A simplified way to test the validity of the SM in the Higgs sector is provided by the so called 'kappa framework' [29,30]. In this framework, the couplings of the Higgs boson to the other SM particles involved at leading order (LO) in perturbation theory for the process under study are dressed with scaling factors . In this simplified approach, based on several assumptions described in Section 10.2 of Ref. [29], production and decay yields are scaled by powers of the corresponding coupling modifier defined as the ratio of the coupling between the particle and the Higgs boson to its SM value. Any significant deviation of a measured from unity would indicate the presence of physics beyond the SM in the tested interaction. In this work, only the coupling modifiers , , , and are considered for single-Higgs interactions (in addition to the modifier that impacts the NLO EW corrections as described in the following). They describe the modifications of the SM Higgs boson coupling to up-type quarks, to down-type quarks, to leptons and to vector bosons ( = or ) respectively. In this parameterisation the interactions between the Higgs boson and the gluons and photons are resolved in terms of the coupling modifiers of the SM particles that enter the loop-level diagrams. New particles contributing to these diagrams are not considered. The total width of the Higgs boson is also parameterised in terms of the coupling modifiers of the individual SM particles, assuming no beyond-the-SM contributions. For double-Higgs production the coupling modifiers , , and 2 are considered. The last of these is related to the interaction vertex, which can be tested in double-Higgs vector-boson fusion (VBF) production (VBF ) as described in the following.
Double-Higgs production is directly sensitive to the Higgs boson self-coupling, starting at the lowest order in perturbation theory. In the SM, the gluon-gluon fusion process (ggF ) accounts for more than 90% of the Higgs boson pair-production → cross-section. The next most abundant process is VBF production, while very small contributions are expected from double-Higgs production in association with a vector boson ( ) and in association with top-quarks (¯). An overview of double-Higgs production at the LHC can be found in Ref. [31].
At lowest order in perturbation theory, the ggF process proceeds via two amplitudes: the first (A 1 ) represented by diagram (a) in Figure 1, and the second (A 2 ) represented by diagram (b). The A 1 amplitude is proportional to the square of the Higgs boson coupling to the top-quark, which scales as 2 , and the A 2 amplitude is proportional to the product of and the Higgs boson self-coupling modifier .
In the SM, the interference between these two amplitudes is destructive and yields an overall cross-section of SM ggF ( → ) = 31.0 +2.1 −7.2 fb at √ = 13 TeV, calculated at NLO in QCD with the measured value of the top-quark mass and corrected to next-to-next-to-leading order (NNLO) including finite top-quark mass effects [30,[32][33][34][35][36][37][38][39][40][41]. The large negative uncertainty originates from the scheme and scale choice of the virtual top-quark mass [41]. Deviations of the ggF cross-section from the SM prediction can therefore be parameterised in terms of the two coupling modifiers and following the prescription described in Refs. [30,[34][35][36][37][38][39][40]. Higher-order QCD corrections do not add further or vertices to the diagrams shown in Figure 1, implying that this parameterisation is applicable to any order in QCD (i.e. also when the amplitudes A 1 and A 2 are modified to include their higher-order QCD corrections). Signal samples for ggF double-Higgs production can be obtained from simulated samples that are generated at different values of these couplings and then combined using morphing techniques, as described in Ref. [27]. Detailed validation studies of this procedure can be found in Ref. [42]. In the SM, the -quark loop contribution to the ggF cross-section is negligible [30,[43][44][45], so its contribution is not included in this analysis. The second most abundant SM double-Higgs process is VBF production, with a predicted SM cross-section of 1.72 ± 0.04 fb at 13 TeV [46][47][48]. At LO in perturbation theory, this process depends on several diagrams that involve the interaction of the Higgs boson with the or vector bosons as shown in Figure 1. The three representative diagrams that enter the total amplitude of the VBF process can be parameterised with different combinations of the , and 2 coupling modifiers [49]. The first diagram, shown in Figure 1(c), is proportional to and , the second, shown in Figure 1(d), to 2 and the last one, shown in Figure 1(e) and related to the quartic interaction vertex , to 2 . The VBF production process can therefore be parameterised using six terms derived from the square of the amplitude described above, which scales as a polynomial of , and 2 . The parameterisation of the signal samples, in terms of yields and kinematic properties, for the double-Higgs VBF process as a function of these coupling modifiers is performed using a set of six independent samples generated for different values of , and 2 . The values of , , and 2 for these six samples were chosen to obtain good statistical precision in the region of parameter space where this analysis is sensitive. The validity of this parameterisation was checked with additional VBF signal samples generated with different values of these coupling modifiers.
The ggF process is sensitive to the sign of relative to the top-quark couplings because of interference between different amplitudes whose leading-order Feynman diagrams are depicted in Figure 1. Similarly, the VBF process provides sensitivity to the relative sign between 2 and .
A complementary approach to study the Higgs boson self-coupling is to use single-Higgs processes, as proposed in Refs.
[20-25]. These processes do not depend on at LO, but the Higgs boson self-coupling contributes to the calculation of the complete NLO EW corrections. In particular, contributes to NLO EW corrections via Higgs boson self-energy loop corrections and via additional diagrams, examples of which are shown in Figure 2. Therefore, an indirect constraint on can be extracted by comparing precise measurements of single-Higgs production and decay yields with the SM predictions corrected for the -dependent NLO EW effects. A framework for a global fit to constrain the Higgs boson self-coupling and the other coupling modifiers was proposed in Refs. [20,21]; the model-dependent assumptions of this parameterisation are described in the same references. In the current work, inclusive production cross-sections, decay branching ratios and differential cross-sections are exploited to increase the sensitivity of the single-Higgs analyses to and . The differential information is encoded through the simplified template cross-section (STXS) framework described in Section III.3 of Ref. [50]. The signal yield in a specific decay channel and STXS bin is then proportional to: signal , where and describe respectively the multiplicative corrections to the expected SM Higgs boson production cross-sections in an STXS bin ( SM, ) and each decay-channel branching ratio (B SM, ) as a function of the values of the Higgs boson self-coupling modifier and the LO-inspired modifiers . The ( × ) coefficients take into account the analysis efficiency times acceptance in each production and decay mode.
The functional dependence of ( , ) and ( , ) on and varies according to the production mode, the decay channel and, more strongly for the and production modes, on the STXS bin. A detailed description of the cross-section and decay-rate dependence on is given in Refs. [51,52]. The STXS information from the VBF, , and production modes is exploited here to constrain and . For the ggF production mode, only the inclusive cross-section dependence on is currently available and it was used in this study, while the STXS bin dependence was not considered.
Conversely, the -modifier can affect the Higgs boson production kinematics and thus modify the analysis efficiency times acceptance in a given STXS bin. This residual dependence was evaluated and found to be negligible for single-Higgs processes, as described in Ref. [51]. Thus the single-Higgs selection acceptances and efficiencies are assumed to be constant as a function of in each STXS bin. A detailed description of the parameterisation of the single-Higgs processes as a function of the coupling modifiers used in this Letter can be found in Ref. [52]. The model under discussion does not allow for any new physics beyond that encoded in the aforementioned and parameters. The dependence of the decay branching ratios and the Higgs boson self-energy on is also taken into account for the double-Higgs analyses when combining them with the single-Higgs results.

Data samples and combined analyses
The results, presented in Sections 5 and 6, are obtained using the full Run 2 dataset collected by the ATLAS experiment [53][54][55] from LHC 13 TeV collisions in the 2015-2018 data-taking period. The integrated luminosity corresponds to 126-139 fb −1 , depending on the trigger selection. A two-level trigger system [56] is used to select events. An extensive software suite [57] is used in the reconstruction and analysis of collision and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.
Each input analysis used in the combination is summarised in Table 1. Details about the individual analyses can be found in the references reported in the same table. Each analysis separates the selected events into different kinematic and topological regions, called categories.

Statistical model and systematic uncertainty correlations
The statistical treatment used in this Letter follows the procedures described in Refs. [65,66]. The results are obtained from a likelihood function ( ì, ì ), where ì represents the vector of the parameters of interest (POI) of the model and ì is a set of nuisance parameters, including the systematic uncertainty contributions and background parameters that are constrained by sidebands or control regions in data. The global likelihood function ( ì, ì ) is obtained as the product of the likelihoods of each input analysis. These are, in turn, products of likelihoods computed in the single analysis categories. The results presented in the following sections are based on the profile-likelihood-ratio test statistic Λ( ì, ì ), and 68% as well as 95% CL intervals are derived in the asymptotic approximation [67]. The CL s approach [68] is only used to derive the cross-section upper limits shown in Section 5.
To derive the expected results, Asimov datasets [67] are produced with all the nuisance parameters set to the values derived from the fit to the data and the parameters of interest fixed to the values corresponding to the hypothesis mentioned in the text.
The basic assumption in performing a statistical combination by using the product of the likelihoods is that the analyses being combined are statistically independent. For this reason the event samples used in the single-Higgs and double-Higgs analyses were checked for overlaps. The overlap among the single-Higgs analyses was checked previously in the combination published in Ref.
[26] and found to be negligible. The event overlap among the three double-Higgs analyses combined for the first time for this result was studied and found to be significantly smaller than 0.1%. These analyses are therefore treated as statistically independent. As a last step, the overlap of event samples between the single-Higgs and double-Higgs analyses, which are combined for the first time in this Letter, was investigated. For most of the categories, this overlap is significantly below the 1% level in either the single-Higgs or the double-Higgs channel, and can therefore be neglected. The only exception is the overlap between the → + − and →¯+ − channels, mainly due to the categories in the → + − analysis, which is found to be at the 4% level in the double-Higgs signal regions. The categories in the → + − channel were removed from the combination used to produce the results presented in the following sections.
A complete discussion of the sources of systematic uncertainty considered in the individual analyses is provided in the publications referenced in Table 1. The correlation model adopted for the systematic uncertainties within the single-Higgs combination is described in detail in Ref. [26].
For this Letter, additional correlations of systematic uncertainties between the double-Higgs analyses and between the single-Higgs and double-Higgs combinations were investigated and implemented as needed. In both cases, systematic uncertainties related to the data-taking conditions, such as those associated with pile-up mis-modelling and the integrated luminosity, are considered to be fully correlated among the input searches. Uncertainties related to physics objects used by multiple searches are treated as correlated where appropriate: experimental uncertainties that are related to the same physics object but determined with different methodologies or implemented with different parameterisations are treated as uncorrelated. Theoretical uncertainties of simulated signal and background processes, such as the single-Higgs and double-Higgs production cross-sections, QCD scale, and proton parton distribution functions are treated as correlated where relevant. The experimental uncertainty of the Higgs boson mass measurement [15] is treated as correlated where relevant. Signal theory uncertainties of the single-Higgs and double-Higgs production modes (e.g., missing higher-order QCD corrections, parton shower, parton distribution functions, etc.) are treated as uncorrelated, while the systematic uncertainties of the decay branching ratios are treated as correlated. For the systematic uncertainties that are constrained significantly in the fit to data, the impact of treating them as correlated or uncorrelated in the combined fit was checked. In general, the impact of these different correlation schemes on the exclusion limits is found to be very small, below the 2% level. Since choosing to treat them as uncorrelated gives slightly larger uncertainties for the parameter of interest, this approach was chosen for the results presented in the following sections.
For the double-Higgs analyses, the most important uncertainties are related to background estimates from data-driven methodologies (derived from data sidebands or control regions) and are therefore not correlated with the single-Higgs analyses. The change of the correlation scheme was found to have a negligible impact on the combined double-Higgs results, except for the theoretical uncertainties of the ggF cross-section, where assuming a correlation loosens the limits on the signal strength by 7% and this is therefore adopted.

Double-Higgs combination results
The double-Higgs boson analyses in the¯¯,¯+ − and¯decay channels referenced in Table 1 are combined in order to place constraints on the production cross-section and the Higgs boson's selfcoupling. First, the value of the signal strength , defined as the ratio of the double-Higgs production cross-section, including only the ggF and VBF processes, to its SM prediction of 32.7 fb [30][31][32][33][34][35][36][37][38][39][40]46] is determined. To produce this result the ratio of the ggF to VBF production cross-sections and the relative kinematic distributions are assumed to be as predicted by the SM, and the other minor production modes are neglected.
This combination yields an observed 95% CL upper limit on of 2.4, with an expected upper limit of 2.9 in the absence of production and 4.0 expected in the SM case. The limits on the signal strength obtained from the individual channels and their combination are shown in Figure 3. The best-fit value obtained from the fit to the data is = −0.7 ± 1.3, which is compatible with the SM prediction of unity, with a -value of 0.2. From the same combination, a 95% CL upper limit on ( → ) of 73 fb is derived (where only ggF and VBF processes are considered), compared with an expected limit of 85 fb assuming no production. When deriving the cross-section limits the theoretical uncertainties on the predicted cross-sections are not included. The cross-section limit as a function of the coupling modifier is shown in Figure 4(a). The signal acceptance of the double-Higgs analyses has a strong dependence on the value of (mainly due to its impact on the distribution), determining the shapes of the exclusion limit curve shown in Figure 4(a).
Constraints on the coupling modifiers are obtained by using the values of the test statistic as a function of in the asymptotic approximation and including the theoretical uncertainty of the cross-section predictions. The parameterisation of NLO EW corrections in the Higgs boson decay and self-energy, as well as in single-Higgs backgrounds, is included when deriving these results, although its impact on the constraints is negligible. With these assumptions, the observed (expected) constraints at 95% CL are −0.6 < < 6.6 (−2.1 < < 7.8). The expected constraint is derived using the SM assumption. More results with different assumptions about the other coupling modifiers are given in Section 6.
The combined double-Higgs channels are also sensitive to the VBF process, and hence to the quartic interaction. The 95% CL observed VBF cross-section upper limit as a function of 2 is shown in Figure 4( Obs. Exp. Observed limit Expected limit ( HH = 0 hypothesis) Expected limit ±1 Expected limit ±2

HH bb + + bb + bbbb
Observed limit (95% CL) Expected limit (95% CL) ( HH = 0 hypothesis) Expected limit ±1 Expected limit ±2 Theory prediction SM prediction bb bb + bbbb Combined (b) Figure 4: Observed and expected 95% CL exclusion limits on the production cross-sections of (a) the combined ggF and VBF processes as a function of and (b) the VBF process as a function of 2 , for the three double-Higgs search channels and their combination. The expected limits assume no production or no VBF production respectively. The red line shows (a) the theory prediction for the combined ggF and VBF cross-section as a function of where all parameters and couplings are set to their SM values except for , and (b) the predicted VBF cross-section as a function of 2 . The bands surrounding the red cross-section lines indicate the theoretical uncertainty on the predicted cross-section. The uncertainty band in (b) is smaller than the width of the plotted line.

Single-and double-Higgs combination results
Following the prescriptions described in Section 2 the double-Higgs and single-Higgs analyses summarised in Table 1 are combined to derive constraints on . Several fits to data are performed with different assumptions about the coupling modifiers to other SM particles.
At first, only possible deviations of from its SM value are considered, assuming that all other Higgs boson interactions proceed as predicted by the SM. The values of twice the negative-logarithm of the profile likelihood ratio (−2 ln Λ) as a function of are shown in Figure 5 for the single-Higgs and double-Higgs analyses, and their combination.   is also superimposed (green curve). The observed best-fit value of for the generic model is shifted slightly relative to the other models because of its correlation with the best-fit values of the , and parameters, which are slightly below, but compatible with unity.
The combined observed (expected) constraints obtained under this hypothesis are −0.4 < < 6.3 (−1.9 < < 7.6) at 95% CL. All the expected constraints reported in this section are derived from an Asimov dataset generated for the SM assumption that corresponds to all coupling modifiers equal to unity. The result is driven by the double-Higgs combination as can be seen in Figure 5. The expected test statistic (−2 ln Λ) curve in Figure 5(b) exhibits a 'two-minima-like' structure due to the quadratic dependence of the observed signal yields on the parameter of interest (partially resolved by the kinematic information used in the fit). The observed curve is more parabolic because the best-fit value of is close to the value where the predicted double-Higgs cross-section, shown in Figure 4(a), reaches its minimum.
The main advantage of adding the single-Higgs analyses is the possibility of relaxing assumptions about modifiers for couplings to other SM particles. First, the assumption about the Higgs boson to top-quark coupling modifier, , can be released. Thanks to the strong constraints on from the single-Higgs measurements, the constraints on obtained from a fit with a floating value of are almost as strong as those obtained with its value fixed to unity, as reported in Table 2. Two-dimensional contours of −2 ln Λ in the -plane are shown in Figure 6. All other coupling modifiers are fixed to unity in this fit.   The most generic model allows all of the coupling modifiers , , , , and implemented in this parameterisation to float freely in the fit. The exception is 2 , which is fixed to unity since there is no complete parameterisation of single-Higgs NLO EW corrections as a function of this coupling modifier. A recent work [69], shows that a consistent parameterisation of the and 2 coupling modifiers seems to be possible, though the sensitivity of single-H processes to 2 is shown to be very small.
In the combination of the single-Higgs and double-Higgs analyses, an observed (expected) exclusion of −1.4 < < 6.1 ( −2.2 < < 7.7) is obtained at 95% CL in this less model-dependent fit. The values of all the other coupling modifiers agree with the SM prediction within uncertainties. The values of the test statistic as a function of for this generic model are also shown in Figure 5. It was checked that for a generic model in which 2 also floats freely in the double-Higgs parameterisation, the observed exclusion constraints on weaken by less than 5%. In this approach, the vertex is parameterised in terms of the 2 coupling modifier for the VBF process but the single-Higgs NLO EW corrections are not.

Conclusion
Single-and double-Higgs boson analyses based on the complete LHC Run 2 dataset of 13 TeV proton-proton collisions collected with the ATLAS detector are combined to investigate the Higgs boson self-interaction and shed more light on the Higgs boson potential, the source of EW symmetry breaking in the SM.
Using the three most sensitive double-Higgs decay channels,¯¯,¯+ − and¯, an observed (expected) upper limit of 2.4 (2.9) at 95% CL is set on the double-Higgs signal strength, defined as the sum of the ggF and VBF production cross-sections normalised to its SM prediction. These processes are directly sensitive to the Higgs boson self-coupling. This combination can also be used to set a constraint of −0.6 < < 6.6 at 95% CL on the Higgs boson self-coupling modifier, assuming that the other Higgs boson interactions are as predicted by the SM.
Using the VBF process, a constraint on the 2 coupling modifier of 0.1 < 2 < 2.0 is also derived at 95% CL, assuming all other Higgs boson interactions are as predicted by the SM.
The measurements from the three double-Higgs decay channels are combined with single-Higgs boson cross-section measurements from the the , * , * , + − and¯decay channels to derive constraints on that are either more stringent or less model-dependent. Using this combination and assuming that is the only source of physics beyond the SM, values of outside the range −0.4 < < 6.3 are excluded at 95% CL, with an expected excluded range of −1.9 < < 7.6. If assumptions about the other coupling modifiers, , , and , are relaxed, this constraint becomes −1.4 < < 6.1 at 95% CL, where the expected interval under the SM assumption is −2.2 < < 7.7. This constraint on the Higgs boson self-coupling is not quite as strong but less model-dependent. This study provides the most stringent constraints on Higgs boson self-interactions to date.