Combination of searches for Higgs boson pairs in pp collisions at root s=13 TeV with the ATLAS detector

This letter presents a combination of searches for Higgs boson pair production using up to 36.1fb − 1 of proton–proton collision data at a centre-of-mass energy √ s = 13 TeV recorded with the ATLAS detector at the LHC. The combination is performed using six analyses searching for Higgs boson pairs decaying into the b ¯ bb ¯ b , b ¯ bW + W − , b ¯ b τ + τ − , W + W − W + W − , b ¯ b γγ and W + W − γγ ﬁnal states. Results are presented for non-resonant and resonant Higgs boson pair production modes. No statistically signiﬁcant excess in data above the Standard Model predictions is found. The combined observed (expected) limit at 95% conﬁdence level on the non-resonant Higgs boson pair production cross-section is 6 . 9 (10) times the predicted Standard Model cross-section. Limits are also set on the ratio ( κ λ ) of the Higgs boson self-coupling to its Standard Model value. This ratio is constrained at 95% conﬁdence level in observation (expectation) to − 5 . 0 < κ λ < 12 . 0 ( − 5 . 8 < κ λ < 12 . 0). In addition, limits are set on the production of narrow scalar resonances and spin-2 Kaluza–Klein Randall–Sundrum gravitons. Exclusion regions are also provided in the parameter space of the habemus Minimal Supersymmetric Standard Model and the Electroweak Singlet Model.


Introduction
The discovery of the Higgs boson (H ) [1,2] at the Large Hadron Collider (LHC) [3] in 2012 has experimentally confirmed the Brout-Englert-Higgs (BEH) mechanism of electroweak symmetry breaking and mass generation [4][5][6]. The BEH mechanism not only predicts the existence of a massive scalar particle, but also requires this scalar particle to couple to itself. Therefore, observing the production of Higgs boson pairs (H H) and measuring the Higgs boson self-coupling λ H H H is a crucial validation of the BEH mechanism.
Any deviation from the Standard Model (SM) predictions would open a window to new physics. Moreover, the form of the Higgs field potential, which generates the Higgs boson self-coupling after electroweak symmetry breaking, can have important cosmological implications, involving, for example, predictions for vacuum stability or models in which the Higgs boson acts as the inflation field [7][8][9][10].
In the SM, the gluon-gluon fusion pp → H H process (ggF) accounts for more than 90% of the Higgs boson pair production cross-section, and only this production mode is considered here. It proceeds via two amplitudes: the first (A 1 ) represented by the diagrams (a) and (b), and the second (A 2 ) represented by the diagram (c) in Fig. 1. The interference between these two amplitudes is de-E-mail address: atlas .publications @cern .ch. structive and yields an overall cross-section of σ SM ggF (pp → H H) = 33.5 +2. 4 −2.8 fb at √ s = 13 TeV [11], calculated first at next-to-leading order (NLO) in QCD with the heavy top-quark approximation [12], then numerically with full top-quark mass dependence [13] (confirmed later in Ref. [14] and analytically computed with some approximation in Ref. [15]) corrected at next-to-next-to-leading order (NNLO) [16] in QCD matched with next-to-next-to-leading logarithmic (NNLL) resummation in the heavy top-quark limit [17,18]. The Higgs boson mass used in these calculations and for all results in this paper is m H = 125.09 GeV [19]. Beyond-the-Standard-  [21,22], constraining κ λ to the range −8 < κ λ < 14 at 95% confidence level (CL). The Higgs boson self-coupling is discussed in the context of BSM models in Refs. [22,23].
Several BSM models also predict the existence of heavy particles decaying into a pair of Higgs   troweak Singlet Models (EWK-singlet) [11,[29][30][31] predict, in addition to the Higgs boson, a second, heavier, CP-even scalar that can decay into two SM Higgs bosons. In the EWK-singlet model, the scalar states are mixed, with a mixing angle α. The ratio of the vacuum expectation value of the additional singlet to that of the SM Higgs doublet, tan β, is a free parameter. In the hMSSM, the CP-even states also mix, and the model's phenomenology can be described by the mass (m A ) of a third, CP-odd, resonance and the ratio of the vacuum expectation values of the two Higgs doublets, tan β. Alternatively, the Higgs boson pair can be produced resonantly through the decay of a spin-2 Kaluza-Klein (KK) graviton, as predicted in the Randall-Sundrum (RS) model of warped extra dimensions [32]. A schematic diagram for production of a heavy resonance followed by its decay into a Higgs boson pair is shown in Fig. 1(d).
This letter presents a combination of results from searches for both non-resonant and resonant Higgs boson pair production in proton-proton (pp) collisions at

Analysis description
The analysis strategies for each of the final states considered in this letter are summarised below.
• The bbbb analysis is performed using four anti-k t jets reconstructed with a radius parameter R = 0.4 [51,52] (resolved analysis) or two large-R jets with R = 1.0 (boosted analysis).
The dataset of the resolved analysis is split according to the years 2015 and 2016, and then statistically combined taking into account the different trigger algorithms used in 2015 and 2016. In part of the 2016 data period, inefficiencies in the online vertex reconstruction affected b-jet triggers that were used in the resolved analysis, reducing the total available integrated luminosity to 27.5 fb −1 . The boosted analysis searches for two large-R jets containing the b-quark pairs from the de-cays of the two Higgs bosons. The large-R jets are identified as originating from b-quarks using a b-tagging algorithm applied to R = 0.2 track-jets [53] associated with the large-R jet [54]. The analysis is divided into three categories: the first category selects events in which each of the two large-R jets has one b-tagged track-jet; the second category requires that one large-R jet contains two b-tagged track-jets and the other large-R jet contains one b-tagged track-jet; the third category requires that both large-R jets contain two b-tagged track-jets. For the SM H H search, only the resolved analysis is used, with two categories, one for the 2015 and another for the 2016 dataset. The resonant search is instead performed with the resolved analysis for masses in the range 260-1400 GeV, with the boosted analysis for masses in the range 800-3000 GeV, and with the combination of the two for masses in the overlapping range 800-1400 GeV.
• The bbW + W − analysis looks for the W W → νqq decay channel, where is an electron or muon, and q is a u, d, s, c quark or anti-quark. The bb pair is selected from two R = 0.4 jets (resolved analysis) or one R = 1.0 large-R jet (boosted analysis), while the jets from the W decay are reconstructed with R = 0.4 jets. The resolved analysis is used in the SM H H search, in the search for a scalar resonance with a mass between 500 and 1400 GeV, and in the search for a KK graviton in the mass range 500 to 800 GeV. The boosted analysis looks for scalar resonances in the mass range 1400 to 3000 GeV and for KK gravitons between 800 and 3000 GeV. The resolved and boosted analyses each use one category. The two analyses are not statistically combined due to a significant overlap between the two signal regions.
• The bbτ + τ − analysis looks for final states with two R = 0.4 b-tagged jets and two τ -leptons. One of the two τ -leptons of the τ + τ − pair is required to decay hadronically, while the other decays either hadronically (τ had τ had ) or leptonically (τ lep τ had ). In the τ lep τ had channel, events are triggered by single lepton triggers (SLT), requiring an electron or a muon in the final state, or by the coincidence of a lepton trigger with a hadronic τ trigger (LTT). In the τ had τ had channel, events are triggered by single hadronic τ triggers (STT) or double hadronic τ triggers (DTT). The analysis is divided into three categories: one selects τ had τ had events, a second selects τ lep τ had events triggered by the SLT, and a third selects τ lep τ had events triggered by the LTT. The τ had τ had and the SLT τ lep τ had categories are used for all model interpretations, while the LTT m H = 125.09 GeV. "L int " indicates the integrated luminosity of the dataset used in the analysis. "Categories" indicates the number of signal categories. "Discriminant" indicates the distribution used in the final limit-setting fit ("c.e." stands for counting events and indicates that a simple event counting was used in the final fit rather than a distribution shape). "Model" indicates which models each analysis tested: NR stands for SM H H signal model, S for a spin-0 scalar model, and G for a KK graviton model. "m S/G " gives the probed mass range for the resonant search. "Ref." reports the reference to the individual final state papers. ν ν 4q, ν ν ν 2q, and ν ν ν ν, with being an electron or muon, q a quark, and ν a neutrino. The q momentum is reconstructed from R = 0.4 jets. In order to suppress Z + jets and tt background, dilepton events are required to have two leptons of the same charge. Events are categorised according to the lepton flavour (ee, eμ and μμ). Three-lepton events are selected if the sum of the lepton charges is ±1. They are divided into two categories according to the number of sameflavour, opposite-charge (SFOS) lepton pairs; one category selects zero SFOS lepton pairs and a second category selects one or two SFOS lepton pairs. Four-lepton events are categorised according to the number of SFOS lepton pairs and the invariant mass (m 4 ) of the four-lepton system. Four categories are defined, requiring that the number of SFOS lepton pairs is less than two or equal to two, and m 4 is smaller or larger than 180 GeV. A total of nine categories are fit simultaneously in the searches for both non-resonant and resonant H H production.
• The bbγ γ analysis searches for a H H pair decaying into bb and γ γ . Two high-p T isolated photons are required to have E T /m γ γ > 0.35 and 0.25 respectively. The events are then analysed using two selections: a 'loose selection' requiring a jet with p T > 40 GeV and a second jet with p T > 25 GeV, and a 'tight selection' where the two jets are required to have p T larger than 100 and 30 GeV. All jets have a radius parameter R = 0.4. Both selections are subdivided into two categories requiring one b-tagged jet or two b-tagged jets. The tight selection is used in the SM H H search and in the search for resonances with masses higher than 500 GeV, while the loose selection is used in the κ λ analysis and in the search for resonances with masses smaller than 500 GeV. The analysis is therefore divided into four categories, but only two of them are simultaneously fit to extract each result.
• The W + W − γ γ analysis searches for a H H pair decaying into γ γ and W W . The analysis uses the same photon selection as the bbγ γ channel and looks for one W decaying leptonically and a second W decaying hadronically (W W → νqq). The hadronic W decay is reconstructed from R = 0.4 jets. Only one category is used in the searches for both non-resonant and resonant H H production.
A summary of the main analysis characteristics is given in Table 1. All analyses impose a series of sequential requirements on kinematic variables to select signal events and suppress back-

Statistical treatment
The statistical interpretation of the combined search results is based on a simultaneous fit to the data for the cross-section of the signal process and nuisance parameters that encode statistical and systematic uncertainties, using the CL S approach [56]. The asymptotic approximation [57] is used in the analysis of all final states and their combination.
All signal regions considered in the simultaneous fit are either orthogonal by construction or have negligible overlap. The overlap due to object misidentification between bbγ γ and bbτ + τ − , and between W + W − γ γ and W + W − W + W − , which are not orthogonal by construction, is evaluated by running the signal and data samples from each channel through the analysis selection of each other channel. Less than 0.1% of simulated signal events overlap between analyses, and no overlap is found in data. There is some irreducible contamination from bbW + W − and bb Z Z events with τ 's in the final state passing the bbτ + τ − selection. This contamination is less than 8% of the bbτ + τ − selected events, and it is not taken into account in the bbτ + τ − analysis, note that including this contribution would increase the analysis sensitivity therefore the results obtained here are slightly conservative. The detector systematic uncertainties, such as those in jet reconstruction, b-jet tagging, electron, muon and photon reconstruction and identification, as well as the uncertainty on the integrated luminosity [58], are correlated across all final states. Uncertainties on the signal acceptance derived by varying the renormalisation and factorisation scales, the parton distribution functions and the parton shower are correlated too. Theoretical and modelling systematic uncertainties of the backgrounds derived using simulated events are not correlated across different analyses because the overlap among their contributions to the different analyses is negligible.
and their statistical combination. The column "Obs." lists the observed limits, "Exp." the expected limits with all statistical and systematic uncertainties, and "Exp. stat." the expected limits obtained including only statistical uncertainties in the fit.

Combination of results on non-resonant Higgs boson pair production
The SM H H analyses use signal samples generated at next-toleading order (NLO) in QCD with Madgraph5_aMC@NLO [59] using the CT10 NLO parton distribution function (PDF) set [60]. Parton showers and hadronisation were simulated with Herwig++ [61] using parameter values from the UE-EE-5-CTEQ6L1 tune [62]. The so-called FTApprox method [63] is applied in the event generation to include finite top-quark mass effects in the real-radiation NLO corrections. The virtual-loop corrections are realised with Higgs effective field theory (HEFT) assuming infinite top-quark mass. The generated events are then corrected with a generator level bin-bybin reweighting of the m H H distribution, which is calculated with finite top-quark mass in full NLO corrections [13]. The branching fractions of the Higgs boson are assumed to be equal to the SM predictions [11]. The limits are determined assuming that all kinematic properties of the H H pair are those predicted by the SM, particularly the m H H distribution, and only the total ggF production cross-section, σ ggF (pp → H H), is allowed to deviate from its SM value. The theoretical uncertainties of σ SM ggF are less than 10% [11] and are not included in the fit results.
The upper limits at 95% CL on the cross-section of the ggF Higgs boson pair production normalised to σ SM ggF are shown in Fig. 2 for the individual final states and their combination. The upper limit for each final state is obtained from a fit with minimal changes from previously published results. The changes include an update of the ggF Higgs boson pair production cross-section from 33.4 fb to 33.5 fb for all final states. Additionally, the bbτ + τ − final state included theoretical uncertainties on the ggF inclusive crosssection, σ SM ggF , which are not considered in the present treatment, and the bbγ γ final state is updated to use an asymptotic approximation to calculate the observed limit instead of the pseudoexperiment method used for its publication. This results in a 10% change in the observed limit of bbγ γ . Moreover, the impact of the asymptotic approximation on all final states combined is found to be 5%.
The combined observed (expected) upper limit on the SM H H production is 6.9 (10) × σ SM ggF (pp → H H , remaining however within the two 2σ uncertainty interval. Detailed comparisons can be found in Ref. [64]. The impact of the systematic uncertainties has been evaluated by recomputing the limit without their inclusion. The limit is then reduced by 13% when removing all of them. The main sources of systematic uncertainty are the modelling of the backgrounds, the statistical uncertainty of simulated events and the τ -lepton reconstruction and identification. When removed the limit reduces by 5%, 3% and 2%, respectively.

Constraints on the Higgs boson self-coupling
The results in Fig. 2 Fig. 3(a), where, with the exception of a small excess in the region below 300 GeV [36] and a small deficit in the 500-600 GeV region, good agreement between data and the as can be seen in Fig. 3. The sensitivity of this analysis is instead affected by a variation in the signal acceptance by up to a factor of three over the probed range of κ λ -values, as shown in Fig. 4(a).
• In the bbγ γ final state, the loose selection is used in the κ λ -scan analysis because the average transverse momentum of the Higgs bosons is lower at large values of κ λ , where |A 2 | 2 dominates the production cross-section. As in the SM H H search, the statistical analysis is performed using the m γ γ distribution, which does not depend on κ λ . The signal acceptance varies by about 30% over the probed range of κ λ -values, as shown in Fig. 4(a). In the previously published analysis [40], LO samples were used for the computation of the signal acceptance, while in this paper the NLO reweighted samples are used, as described above. The signal acceptance times efficiency as a function of κ λ is shown in Fig. 4(a). Given that, for each final state, the same selection was applied over the full scanned κ λ range, the shape of the acceptance times efficiency curve is determined by the variation of the event kinematics as a function of κ λ . For high values of |κ λ | the A 2 term dominates the total amplitude, causing a softer m HH spectrum, and thus a lower acceptance times efficiency. Around κ λ = 2.4 the interference between A 1 and A 2 amplitudes is maximal, producing the hardest m HH spectrum and, consequently, the highest signal acceptance times efficiency.
In Fig. 4(b) the shape of the upper-limit curves approximately follows the inverse of the signal acceptance shown in Fig. 4(a). In the bbbb analysis, the observed limits are more stringent than the values, and thus the excess in data below 300 GeV leads to the observed limits being less stringent than expected. In the bbτ + τ − final state the observed limits are more stringent than the expected limits over the whole range of κ λ , due to a deficit of data relative to the background predictions at high values of the BDT score. The bbγ γ limit shows a weaker dependence on κ λ than the bbbb and bbτ + τ − limits because the bbγ γ acceptance varies less as function of κ λ .
The 95% CL allowed κ λ intervals are given in Table 2. The systematic uncertainties weaken the κ λ limits by less than 10% relative to those obtained with only statistical uncertainties. The final state least (most) affected by systematic uncertainties is bbγ γ (bbbb). The Higgs boson branching fraction depends on κ λ due to Table 2 Allowed κ λ intervals at 95% CL for the bbbb, bbτ + τ − and bbγ γ final states and their combination. The column "Obs." lists the observed results, "Exp." the expected results obtained including all statistical and systematic uncertainties in the fit, and "Exp. stat." the expected results obtained including only the statistical uncertainties.
The effect of non-SM Higgs decay branching fractions due to κ λ variations is not taken into account, which impacts the κ λ intervals by no more than 7%. . This dependence is neglected in the present treatment, but its overall impact on the allowed κ λ interval is evaluated to be no more than 7%. Theory uncertainties on the signal cross section shown in Fig. 4(b) are not taken into account when computing the κ λ limits in Table 2, they affect the limit by less than 8%.

Combination of results for resonant Higgs boson pair production
The resonance decaying into a pair of Higgs bosons is assumed to be either a heavy spin-0 scalar particle, S, with a narrow width or a spin-2 KK graviton, G KK .
The search for the heavy scalar particle S is performed with all six final states included in this combination. With the exception of bbτ + τ − and bbbb, all signal samples were simulated at NLO with MadGraph5_aMC@NLO using the CT10 PDF set. The matrix-element generator was interfaced to Herwig++ with the UE-EE-5-CTEQ6L1 tune. The bbτ + τ − final state uses an LO model generated with MadGraph5_aMC@NLO using the NNPDF 2.3 LO PDF set interfaced to Pythia 8.2 with the A14 tune, while the bbbb final state uses the same LO event generator but interfaced to Her-wig++ with the UE-EE-5-CTEQ6L1 tune.
The scalar resonance search is performed in the mass range 260-3000 GeV, and within this range no statistically significant excess is observed. In the combination, the largest observed deviation from the background expectation is 1σ for the search mass range. The combined upper limit on the cross-section is shown as a func- tion of the resonance mass in Fig. 5(a). Systematic uncertainties have a sizeable effect on the upper limits depending on the probed resonance mass. The total impact of systematics or the impact of a single systematic uncertainty has been evaluated by computing the percentage reduction of the upper limit obtained by removing all systematic uncertainties or a particular source. Overall the systematic uncertainties affect the limit by 12% (11%) for a resonance mass of 1 (3) TeV. Among them, the largest systematic uncertainties are due to the modelling of the backgrounds, impacting the upper limit by 7% (9%) at 1 (3) TeV. The second leading systematic uncertainty comes from b-tagging, that affects the upper limit by 2% at 1 TeV, but its impact is negligible at 3 TeV where relative background and statistical uncertainties increase significantly. At 3 TeV the second leading systematic uncertainty is related to the jet energy scale and resolution, changing the limit by 2%. Interpretations in specific spin-0 BSM models are provided in Section 7.
The search for a spin-2 KK graviton is performed with the bbbb, uncertainties come from b-tagging at low G KK mass, that affect the limit by 3%, and from jet energy scale and resolution at high mass, that affect the upper limit by 2% (3%) at 1 TeV (3 TeV). For  [74] are shown with horizontal lines. The dotted lines indicate the separation between regions where the resonance width is larger than 2% and 5% of the resonance mass. The hatch-marked area corresponds to regions that cannot be excluded because the width of the resonance exceeds 10% of the resonance mass, corresponding to the maximum of the experimental mass resolution among all analysed final states.
For the EWK-singlet model, the experimental limits on the spin-0 resonance (as reported in Section 6) are interpreted as constraints in the m S -sin α plane (where m S is the resonance mass) for tan β = 1 and tan β = 2, shown in Fig. 6(a) and Fig. 6(b) respectively. The expected cross-section for each point in the parameter space is obtained by scaling the heavy Higgs cross-section calculated at NNLO+NNLL [11] with singlet coupling modifiers. The branching fractions are computed with sHDECAY [70]. In this model, the width of the heavy scalar can be large in some regions of the parameter space. Due to the use of narrow-width signal models in the event generation, results presented here are valid only in regions of the model parameter-space where the resonance width ( S ) is smaller than the experimental resolution at the resonance mass. This holds when S /m S < 2% for bbγ γ , S /m S < 5% for bbbb and S /m S < 10% for bbτ + τ − . Therefore, the excluded region in the plot is obtained by combining the three final states for S /m S < 2%, by combining the bbbb and bbτ + τ − final states for 2% < S /m S < 5%, and using only bbτ + τ − for 5% < S /m S < 10%. The hatched region shows points where S /m S ≥ 10%, where no exclusion can be provided. Fig. 7(a) shows limits for the (sin α, tan β) parameter space for m S = 260 GeV where, due to the limited decay phase space, the resonance width is narrow in a wide region of the parameter space.
The experimental limits on a spin-0 resonance are also interpreted as constraints in the m A -tan β plane of the hMSSM model in Fig. 7(b). The expected cross-section for each point in the parameter space is obtained using the gluon-gluon fusion crosssection from SUSHI 1.5.0 [71,72] and the branching fractions computed with HDECAY 6.4.2 [73].
The excluded region is more than doubled along tan β relative to the previous combined results in Ref.
[42] at 8 TeV, and excludes values of m A from 190 GeV to 560 GeV depending on tan β. The kink at low tan β and high m A values is caused by removing the bbγ γ final state from the combination in the region where the predicted width of the heavy CP-even Higgs boson is larger than the experimental resolution on m S in the bbγ γ analysis.

Conclusion
A statistical combination of six final states bbbb, bbW

Acknowledgements
We thank CERN for the very successful operation of the LHC, as well as the support staff from our institutions without whom ATLAS could not be operated efficiently.
We acknowledge the support of ANPCyT,