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

TL;DR

This paper combines multiple searches for Higgs boson pair production at 13 TeV with the ATLAS detector, setting limits on production cross-sections, Higgs self-coupling, and new physics models, with no significant excess observed.

Contribution

It provides the first combined analysis of various Higgs pair decay channels at 13 TeV, improving constraints on Higgs self-coupling and new physics scenarios.

Findings

No significant excess above Standard Model predictions.

Limits set on Higgs pair production cross-section and self-coupling.

Constraints on new physics models like resonances and extra dimensions.

Abstract

This letter presents a combination of searches for Higgs boson pair production using up to 36.1 fb$^{-1}$ of proton-proton collision data at a centre-of-mass energy $\sqrt{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 bbbb, bbWW, bb$\tau\tau$, WWWW, bb$\gamma \gamma$ and WW$\gamma\gamma$ final states. Results are presented for non-resonant and resonant Higgs boson pair production modes. No statistically significant excess in data above the Standard Model predictions is found. The combined observed (expected) limit at 95% confidence 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 ($ \kappa_{\lambda} $) of the Higgs boson self-coupling to its Standard Model value. This ratio is constrained at 95% confidence level in observation (expectation) to $ -5.0 < \kappa_{\lambda} <12.0 $ ($ -5.8 < \kappa_{\lambda} <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.

Figures

Tables

Table 1. Table 1: Summary of the main characteristics of the analyses used in the non-resonant and resonant searches. The resonant analysis characteristics are indicated between square brackets. “ ℬ ​ ( H ​ H → x ​ x ¯ ​ y ​ y ¯ ) ℬ → 𝐻 𝐻 𝑥 ¯ 𝑥 𝑦 ¯ 𝑦 \mathcal{B}(HH\to x\bar{x}y\bar{y}) ” indicates the branching fraction of the H ​ H 𝐻 𝐻 HH pair where one H 𝐻 H decays into x ​ x ¯ 𝑥 ¯ 𝑥 x\bar{x} and the other decays into y ​ y ¯ 𝑦 ¯ 𝑦 y\bar{y} . The branching fraction values are taken from Ref.  [ 11 ] for a Higgs boson mass m H = 125.09 subscript 𝑚 𝐻 125.09 m_{H}=125.09 GeV. “ ℒ i ​ n ​ t subscript ℒ i 𝑛 𝑡 \mathcal{L}_{\mathrm{i}nt} ” 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 𝐻 𝐻 HH signal model, S 𝑆 S for a spin-0 scalar model, and G 𝐺 G for a KK graviton model. “ m S / G subscript 𝑚 𝑆 𝐺 m_{S/G} ” gives the probed mass range for the resonant search. “Ref.” reports the reference to the individual final state papers.

	$b\overline{b}b\overline{b}$	$b\overline{b}{W}^{+}{W}^{-}$	$b\overline{b}{\tau }^{+}{\tau }^{-}$	$W^{+}{W}^{-}{W}^{+}{W}^{-}$	$b\overline{b}\gamma \gamma$	$W^{+}{W}^{-}\gamma \gamma$
$ℬ\left(HH\to x\overline{x}y\overline{y}\right)$	0.34	0.25	0.073	0.046	$2.6\cdot {10}^{-3}$	$1.0\cdot {10}^{-3}$
${ℒ}_{int}$ [fb-1]	27.5 [36.1]	36.1	36.1	36.1	36.1	36.1
Categories	2 [2–5]	1 [1]	3 [2–3]	9 [9]	2 [2]	1 [1]
Discriminant	$m_{HH}$ [$m_{HH}$]	c.e. [$m_{HH}$]	BDT [BDT]	c.e. [c.e.]	$m_{\gamma \gamma }$ [$m_{HH}$]	$m_{\gamma \gamma }$ [$m_{\gamma \gamma }$]
Model	NR [$S/G$]	NR [$S/G$]	NR [$S/G$]	NR [$S$]	NR [$S$]	NR [$S$]
$m_{S/G}$ [TeV]	[0.26–3.00]	[0.50–3.00]	[0.26–1.00]	[0.26–0.50]	[0.26–1.00]	[0.26–0.50]
Ref.	[36]	[37]	[38]	[39]	[40]	[41]

Table 2. Table 2: Allowed κ λ subscript 𝜅 𝜆 \kappa_{\lambda} intervals at 95% CL for the b ​ b ¯ ​ b ​ b ¯ 𝑏 ¯ 𝑏 𝑏 ¯ 𝑏 b\bar{b}b\bar{b} , b ​ b ¯ ​ τ + ​ τ − 𝑏 ¯ 𝑏 superscript 𝜏 superscript 𝜏 b\bar{b}\tau^{+}\tau^{-} and b ​ b ¯ ​ γ ​ γ 𝑏 ¯ 𝑏 𝛾 𝛾 b\bar{b}\gamma\gamma 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 κ λ subscript 𝜅 𝜆 \kappa_{\lambda} variations is not taken into account, which impacts the κ λ subscript 𝜅 𝜆 \kappa_{\lambda} intervals by no more than 7%.

	Allowed ${\kappa }_{\lambda }$ interval at 95% CL
Final state	Obs.			Exp.			Exp. stat.
$b\overline{b}b\overline{b}$	$-10.8805$	—	$20.0808$	$-11.6338$	—	$18.7727$	$-9.79596$	—	$16.2699$
$b\overline{b}{\tau }^{+}{\tau }^{-}$	$-7.40278$	—	$15.739$	$-8.85064$	—	$16.761$	$-7.83268$	—	$15.5231$
$b\overline{b}\gamma \gamma$	$-8.06538$	—	$13.0955$	$-8.09256$	—	$13.1123$	$-7.90249$	—	$12.9228$
Combination	$-4.97183$	—	$12.0122$	$-5.81589$	—	$12.0083$	$-5.28363$	—	$11.507$

Full text

\AtlasTitle

Combination of searches for Higgs boson pairs in $pp$ collisions at $\sqrt{s}=$13 TeV with the ATLAS detector \AtlasAbstract This letter presents a combination of searches for Higgs boson pair production using up to 36.1 ${\mathrm{fb}^{-1}}$ of proton–proton collision data at a centre-of-mass energy $\sqrt{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\bar{b}b\bar{b}$, $b\bar{b}W^{+}W^{-}$, $b\bar{b}\tau^{+}\tau^{-}$, $W^{+}W^{-}W^{+}W^{-}$, $b\bar{b}\gamma\gamma$ and $W^{+}W^{-}\gamma\gamma$ final states. Results are presented for non-resonant and resonant Higgs boson pair production modes. No statistically significant excess in data above the Standard Model predictions is found. The combined observed (expected) limit at 95% confidence level on the non-resonant Higgs boson pair production cross-section is $6.91438$ ($10.2769$) times the predicted Standard Model cross-section. Limits are also set on the ratio ($\kappa_{\lambda}$) of the Higgs boson self-coupling to its Standard Model value. This ratio is constrained at 95% confidence level in observation (expectation) to -4.97183$<\kappa_{\lambda}<$12.0122 (-5.81589$<\kappa_{\lambda}<$12.0083). 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.

\AtlasRefCodeHDBS-2018-58 \AtlasJournalPhys. Lett. B \PreprintIdNumberCERN-EP-2019-099 \AtlasCoverSupportingNoteThe combination and interpretation of $HH$ searches in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detectorhttps://cds.cern.ch/record/2653800/ \AtlasCoverCommentsDeadline \AtlasCoverAnalysisTeam John Alison, Nansi Andari, Florian Beisiegel, Alessandra Betti, Petar Bokan, Elizabeth Brost, Patrick Bryant, Tyler Burch, Leo Cerda, David Delgove, Biagio Di Micco, Jochen Dingfelder, David Englert, Yaquan Fang, Arnaud Ferrari, Douglas Gingrich, Kathryn Grimm, Carl Gwilliam, Jonathan Hays, Tatjana Lenz, Qi Li, Gabriel Palacino, James Robinson, Nikolaos Rompotis, Eleonora Rossi, Jana Schaarschmidt, Chase Shimmin, Magdalena Slawinska, Xiaohu Sun, Alan Taylor, Tulin Varol, David Wardrope, Yu Zhang

\AtlasCoverEdBoardMemberKathrin Becker, John Hobbs, Jianming Qian (chair) \AtlasCoverEgroupEditorsatlas-hdbs-2018-58-editors@cern.ch \AtlasCoverEgroupEdBoardatlas-hdbs-2018-58-editorial-board@cern.ch \AtlasJournalRefPhys. Lett. B 800 (2020) 135103 \AtlasDOI10.1016/j.physletb.2019.135103

1 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 ($HH$) and measuring the Higgs boson self-coupling $\lambda_{HHH}$ 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\to HH$ 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 ($\mathcal{A}_{1}$) represented by the diagrams (a) and (b), and the second ($\mathcal{A}_{2}$) represented by the diagram (c) in Figure 1.

The interference between these two amplitudes is destructive and yields an overall cross-section of $\sigma^{\mathrm{S}M}_{\mathrm{g}gF}(pp\to HH)=33.5^{+2.4}_{-2.8}$ fb at $\sqrt{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-Model (BSM) scenarios can bring substantial enhancement of this cross-section by modifying the relative sign of $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$, and by increasing $\mathcal{A}_{2}$. The $\mathcal{A}_{2}$ amplitude is proportional to the Higgs self-coupling $\lambda_{HHH}$. The Higgs boson self-coupling modifier due to BSM scenarios is defined as $\kappa_{\lambda}=\lambda_{HHH}/\lambda^{\mathrm{S}M}_{HHH}$. In this analysis, all other Higgs boson couplings are assumed to have SM values. Indirect limits on $\kappa_{\lambda}$ have been obtained using the measurements of single Higgs boson production and decay [20] and electroweak precision observables [21, 22], constraining $\kappa_{\lambda}$ to the range $-8<\kappa_{\lambda}<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 bosons. Two-Higgs-Doublet Models [24], models inspired by the Minimal Supersymmetric Standard Model (MSSM) like habemus MSSM (hMSSM)  [25, 26, 27, 28], and Electroweak Singlet Models (EWK-singlet) [29, 30, 31, 11] 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 $\alpha$. The ratio of the vacuum expectation value of the additional singlet to that of the SM Higgs doublet, $\tan\beta$, 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\beta$. 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 Figure 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 $\sqrt{s}=13$ TeV. The data were collected with the ATLAS detector [33, 34, 35] and correspond to an integrated luminosity of up to 36.1 ${\mathrm{fb}^{-1}}$. The combination includes all published ATLAS $HH$ search analyses using $\sqrt{s}=13$ TeV data, namely those studying the final states $b\bar{b}b\bar{b}$ [36], $b\bar{b}W^{+}W^{-}$ [37], $b\bar{b}\tau^{+}\tau^{-}$ [38], $W^{+}W^{-}W^{+}W^{-}$ [39], $b\bar{b}\gamma\gamma$ [40] and $W^{+}W^{-}\gamma\gamma$ [41].

Previous combinations of searches for $HH$ pair production were performed at $\sqrt{s}=8$ TeV by the ATLAS experiment [42] and at $\sqrt{s}=13$ TeV by the CMS experiment [43] combining the final states $b\bar{b}b\bar{b}$ [44, 45, 46, 47], $bb$VV [48], $b\bar{b}\tau^{+}\tau^{-}$ [49] and $b\bar{b}\gamma\gamma$ [50].

2 Analysis description

The analysis strategies for each of the final states considered in this letter are summarised below.

• •

The $b\bar{b}b\bar{b}$ 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 decays 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 $HH$ 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 $b\bar{b}W^{+}W^{-}$ analysis looks for the $WW\to\ell\nu qq$ decay channel, where $\ell$ is an electron or muon, and $q$ is a $u,d,s,c$ quark or anti-quark. The $b\bar{b}$ 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 $HH$ 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 $b\bar{b}\tau^{+}\tau^{-}$ analysis looks for final states with two $R=0.4$ $b$-tagged jets and two $\tau$-leptons. One of the two $\tau$-leptons of the $\tau^{+}\tau^{-}$ pair is required to decay hadronically, while the other decays either hadronically ($\tau_{\mathrm{h}ad}\tau_{\mathrm{h}ad}$) or leptonically ($\tau_{\mathrm{l}ep}\tau_{\mathrm{h}ad}$). In the $\tau_{\mathrm{l}ep}\tau_{\mathrm{h}ad}$ 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 $\tau$ trigger (LTT). In the $\tau_{\mathrm{h}ad}\tau_{\mathrm{h}ad}$ channel, events are triggered by single hadronic $\tau$ triggers (STT) or double hadronic $\tau$ triggers (DTT). The analysis is divided into three categories: one selects $\tau_{\mathrm{h}ad}\tau_{\mathrm{h}ad}$ events, a second selects $\tau_{\mathrm{l}ep}\tau_{\mathrm{h}ad}$ events triggered by the SLT, and a third selects $\tau_{\mathrm{l}ep}\tau_{\mathrm{h}ad}$ events triggered by the LTT. The $\tau_{\mathrm{h}ad}\tau_{\mathrm{h}ad}$ and the SLT $\tau_{\mathrm{l}ep}\tau_{\mathrm{h}ad}$ categories are used for all model interpretations, while the LTT $\tau_{\mathrm{l}ep}\tau_{\mathrm{h}ad}$ category is used in the SM $HH$ search (excluding the $\kappa_{\lambda}$ analysis) and in resonant searches up to a mass of 800 GeV.

• •

The $W^{+}W^{-}W^{+}W^{-}$ analysis looks for channels with leptonic and/or hadronic $W$ decays. Three channels are identified: $\ell\nu\,\ell\nu\,4q$, $\ell\nu\,\ell\nu\,\ell\nu\,2q$, and $\ell\nu\,\ell\nu\,\ell\nu\,\ell\nu$, with $\ell$ being an electron or muon, $q$ a quark, and $\nu$ a neutrino. The $q$ momentum is reconstructed from $R=0.4$ jets. In order to suppress $Z+$jets and $t\bar{t}$ background, dilepton events are required to have two leptons of the same charge. Events are categorised according to the lepton flavour ($ee$, $e\mu$ and $\mu\mu$). Three-lepton events are selected if the sum of the lepton charges is $\pm$1. They are divided into two categories according to the number of same-flavour, 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\ell}$) 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\ell}$ is smaller or larger than 180 GeV. A total of nine categories are fit simultaneously in the searches for both non-resonant and resonant $HH$ production.

• •

The $b\bar{b}\gamma\gamma$ analysis searches for a $HH$ pair decaying into $b\bar{b}$ and $\gamma\gamma$. Two high-$p_{\mathrm{T}}$ isolated photons are required to have $E_{\text{T}}/m_{\gamma\gamma}>$ 0.35 and 0.25 respectively. The events are then analysed using two selections: a ‘loose selection’ requiring a jet with $p_{\mathrm{T}}>40$ GeV and a second jet with $p_{\mathrm{T}}>25$ GeV, and a ‘tight selection’ where the two jets are required to have $p_{\mathrm{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 $HH$ search and in the search for resonances with masses higher than 500 GeV, while the loose selection is used in the $\kappa_{\lambda}$ 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^{-}\gamma\gamma$ analysis searches for a $HH$ pair decaying into $\gamma\gamma$ and $WW$. The analysis uses the same photon selection as the $b\bar{b}\gamma\gamma$ channel and looks for one $W$ decaying leptonically and a second $W$ decaying hadronically ($WW\to\ell\nu 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 $HH$ 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 backgrounds. The $b\bar{b}\tau^{+}\tau^{-}$ analysis uses a boosted decision tree (BDT) [55] distribution as the final discriminant for both the non-resonant and resonant searches. For the resonant searches, the $b\bar{b}b\bar{b}$, $b\bar{b}W^{+}W^{-}$ and $b\bar{b}\gamma\gamma$ analyses use the $HH$ invariant mass ($m_{HH}$) as the final discriminant, the $W^{+}W^{-}\gamma\gamma$ analysis uses the $\gamma\gamma$ invariant mass ($m_{\gamma\gamma}$), while $W^{+}W^{-}W^{+}W^{-}$ uses simple event counting. For the SM $HH$ search, the $b\bar{b}b\bar{b}$ analysis uses the $m_{HH}$ distribution as a discriminant, profiting from the difference between the shapes of the signal and the dominant multi-jet background. The $b\bar{b}\gamma\gamma$ and $W^{+}W^{-}\gamma\gamma$ analyses fit the $m_{\gamma\gamma}$ distribution to extract both the signal yield and the background expectation, while the $b\bar{b}W^{+}W^{-}$ and $W^{+}W^{-}W^{+}W^{-}$ analyses use event counting.

3 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 CLS 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 $b\bar{b}\gamma\gamma$ and $b\bar{b}\tau^{+}\tau^{-}$, and between $W^{+}W^{-}\gamma\gamma$ 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 $b\bar{b}W^{+}W^{-}$ and $b\bar{b}ZZ$ events with $\tau$’s in the final state passing the $b\bar{b}\tau^{+}\tau^{-}$ selection. This contamination is less than 8% of the $b\bar{b}\tau^{+}\tau^{-}$ selected events, and it is not taken into account in the $b\bar{b}\tau^{+}\tau^{-}$ 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.

4 Combination of results on non-resonant Higgs boson pair production

The SM $HH$ analyses use signal samples generated at next-to-leading 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-by-bin reweighting of the $m_{HH}$ 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]. For the SM $HH$ search, upper limits are extracted for the cross-section $\sigma_{\mathrm{g}gF}(pp\to HH)$ of $HH$ production and are normalised by the SM $HH$ cross-section $\sigma^{\mathrm{S}M}_{\mathrm{g}gF}(pp\to HH)$. The limits are determined assuming that all kinematic properties of the $HH$ pair are those predicted by the SM, particularly the $m_{HH}$ distribution, and only the total ggF production cross-section, $\sigma_{\mathrm{g}gF}(pp\to HH)$, is allowed to deviate from its SM value. The theoretical uncertainties of $\sigma_{\mathrm{g}gF}^{\mathrm{S}M}$ 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 $\sigma_{\mathrm{g}gF}^{\mathrm{S}M}$ are shown in Figure 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 $b\bar{b}\tau^{+}\tau^{-}$ final state included theoretical uncertainties on the ggF inclusive cross-section, $\sigma_{\mathrm{g}gF}^{\mathrm{S}M}$, which are not considered in the present treatment, and the $b\bar{b}\gamma\gamma$ final state is updated to use an asymptotic approximation to calculate the observed limit instead of the pseudo-experiment method used for its publication. This results in a 10% change in the observed limit of $b\bar{b}\gamma\gamma$. 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 $HH$ production is $6.91438$ ($10.2769$) $\times$ $\sigma^{\mathrm{S}M}_{\mathrm{g}gF}(pp\to HH)$. The expected limit is similar to the CMS result of 12.8 $\times$ $\sigma^{\mathrm{S}M}_{\mathrm{g}gF}(pp\to HH)$. The observed limit is more stringent for the ATLAS result than the CMS result of 22.2 $\times$ $\sigma^{\mathrm{S}M}_{\mathrm{g}gF}(pp\to HH)$ because the three leading channels ($b\bar{b}b\bar{b}$, $b\bar{b}\tau^{+}\tau^{-}$ and $b\bar{b}\gamma\gamma$) have a data deficit in ATLAS and an excess in CMS [43], remaining however within the two 2$\sigma$ 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 $\tau$-lepton reconstruction and identification. When removed the limit reduces by 5%, 3% and 2%, respectively.

5 Constraints on the Higgs boson self-coupling

The results in Figure 2 show that the sensitivity of the SM $HH$ search is driven by the final states $b\bar{b}b\bar{b}$, $b\bar{b}\tau^{+}\tau^{-}$ and $b\bar{b}\gamma\gamma$. These final states are used to set constraints on the Higgs boson self-coupling modifier $\kappa_{\lambda}=\lambda_{HHH}/\lambda_{HHH}^{\mathrm{S}M}$. After setting all couplings to fermions and bosons to their SM values, a scan of the self-coupling modifier $\kappa_{\lambda}$ is performed. The $\kappa_{\lambda}$ factor affects both the production cross-section and the kinematic distributions of the Higgs boson pairs, by modifying the $\mathcal{A}_{2}$ production amplitude. It can also affect the Higgs boson branching fractions due to NLO electroweak corrections [20], but this dependence is neglected in the following.

The signal used in the $\kappa_{\lambda}$ fit was simulated according to the following procedure. For each value of $\kappa_{\lambda}$ the $m_{HH}$ spectrum is computed at the generator-level, using the leading-order (LO) version of MadGraph5_aMC@NLO [59] with the NNPDF 2.3 LO [65] PDF set, together with Pythia 8.2 [66] for the showering model using the A14 tune [67]. Because only one amplitude of Higgs boson pair production depends on $\kappa_{\lambda}$, linear combinations of three LO samples generated with different values of $\kappa_{\lambda}$ are sufficient to make predictions for any value of $\kappa_{\lambda}$. Binned ratios of the $m_{HH}$ distributions to the SM distribution are computed for all $\kappa_{\lambda}$ values and then used to reweight the events of NLO SM $HH$ signal samples, generated using the full detector simulation. This procedure is validated by comparing kinematic distributions obtained with the reweighting procedure applied to the LO SM sample and LO samples generated with the actual $\kappa_{\lambda}$ values set in the event generator. The two sets of distributions are found to be in agreement. This procedure assumes that higher order QCD corrections on the differential cross-section as a function of $m_{HH}$ are independent of $\kappa_{\lambda}$. The reweighted NLO signal sample is used to compute the signal acceptance and the kinematic distributions for different values of $\kappa_{\lambda}$.

This letter presents $\kappa_{\lambda}$ results for the first time in the ATLAS $b\bar{b}b\bar{b}$ and $b\bar{b}\tau^{+}\tau^{-}$ final states and incorporates the previously published result for the $b\bar{b}\gamma\gamma$ final state. The $\kappa_{\lambda}$ analyses closely follow the SM $HH$ search, with some exceptions which are discussed below for each final state.

• •

In the $b\bar{b}b\bar{b}$ final state, the same analysis selection and final discriminant are used in the $\kappa_{\lambda}$-scan analysis and in the SM $HH$ search. The distribution of the final discriminant $m_{HH}$ is shown in Figure 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 GeVregion, good agreement between data and the expected background is observed. The shape of the $m_{HH}$ distribution has a strong dependence on $\kappa_{\lambda}$, and the signal acceptance varies by a factor 2.5 over the probed range of $\kappa_{\lambda}$-values ($-20\leq\kappa_{\lambda}\leq 20$) shown in Figure 4(a). The two effects together determine how the exclusion limits on the cross-section of the $HH$ production vary as a function of $\kappa_{\lambda}$.

• •

In the $b\bar{b}\tau^{+}\tau^{-}$ final state, as in the SM $HH$ search, both $\tau_{\mathrm{l}ep}\tau_{\mathrm{h}ad}$ and $\tau_{\mathrm{h}ad}\tau_{\mathrm{h}ad}$ events are used. In contrast with the SM $HH$ search, LTT $\tau_{\mathrm{l}ep}\tau_{\mathrm{h}ad}$ events (see Section 2) are not used given their negligible contribution. The SM $HH$ search and the $\kappa_{\lambda}$-scan analysis use the same sets of variables to build BDT discriminants. For the $\kappa_{\lambda}$-scan the BDTs are retrained using the NLO SM signal sample reweighted with $\kappa_{\lambda}=20$, ensuring good sensitivity over the whole range of probed $\kappa_{\lambda}$-values. The BDT score distributions are used in the fit to compute the final results. The shape of the $b\bar{b}\tau^{+}\tau^{-}$ BDT distributions does not show a $\kappa_{\lambda}$ dependence as strong as in the $b\bar{b}b\bar{b}$ final state, as can be seen in Figure 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 $\kappa_{\lambda}$-values, as shown in Figure 4(a).

• •

In the $b\bar{b}\gamma\gamma$ final state, the loose selection is used in the $\kappa_{\lambda}$-scan analysis because the average transverse momentum of the Higgs bosons is lower at large values of $\kappa_{\lambda}$, where $|\mathcal{A}_{2}|^{2}$ dominates the production cross-section. As in the SM $HH$ search, the statistical analysis is performed using the $m_{\gamma\gamma}$ distribution, which does not depend on $\kappa_{\lambda}$. The signal acceptance varies by about 30% over the probed range of $\kappa_{\lambda}$-values, as shown in Figure 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.

\floatsetup

[figure]style=plain,subcapbesideposition=center

The signal acceptance times efficiency as a function of $\kappa_{\lambda}$ is shown in Figure 4(a). Given that, for each final state, the same selection was applied over the full scanned $\kappa_{\lambda}$ range, the shape of the acceptance times efficiency curve is determined by the variation of the event kinematics as a function of $\kappa_{\lambda}$. For high values of $|\kappa_{\lambda}|$ the $\mathcal{A}_{2}$ term dominates the total amplitude, causing a softer $m_{\mathrm{H}H}$ spectrum, and thus a lower acceptance times efficiency. Around $\kappa_{\lambda}$ = 2.4 the interference between $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ amplitudes is maximal, producing the hardest $m_{\mathrm{H}H}$ spectrum and, consequently, the highest signal acceptance times efficiency.

In each analysis, and in their combination, the 95% CL upper limit on the $\sigma_{\mathrm{g}gF}(pp\to HH)$ cross-section is computed for different values of $\kappa_{\lambda}$. The results are shown in Figure 4(b). The theoretical $\sigma_{\mathrm{g}gF}(pp\to HH)$ cross-section as a function of $\kappa_{\lambda}$ is overlaid in the figure. It is computed by multiplying the $HH$ SM cross-section $\sigma^{\mathrm{S}M}_{\mathrm{g}gF}(pp\to HH)$ by the ratio $R\left(\kappa_{\lambda}\right)$ of the $pp\to HH$ cross-section computed at $\kappa_{\lambda}$, $\sigma_{\mathrm{g}gF}^{\kappa_{\lambda}}(pp\to HH)$ to the same quantity $\sigma_{\mathrm{g}gF}^{\kappa_{\lambda}=1}(pp\to HH)$ computed at $\kappa_{\lambda}=1$. The $R(\kappa_{\lambda})$ factor is computed at NNLO+NNLL with the infinite top-quark mass approximation  [68]. The resulting observed (expected) confidence interval at 95% CL for $\kappa_{\lambda}$ is: $-4.97183$ $<\kappa_{\lambda}<$ $12.0122$ ($-5.81589$ $<\kappa_{\lambda}<$ $12.0083$) .

In Figure 4(b) the shape of the upper-limit curves approximately follows the inverse of the signal acceptance shown in Figure 4(a). In the $b\bar{b}b\bar{b}$ analysis, the observed limits are more stringent than the expected limits at low values of $\kappa_{\lambda}$. For these $\kappa_{\lambda}$ values the signal $m_{HH}$ distributions have significant populations in the region 500-600 GeV, where the data deficit sits, as explained above. For larger values of $\kappa_{\lambda}$ the $m_{HH}$ distribution is shifted to lower $m_{HH}$ values, and thus the excess in data below 300 GeV leads to the observed limits being less stringent than expected. In the $b\bar{b}\tau^{+}\tau^{-}$ final state the observed limits are more stringent than the expected limits over the whole range of $\kappa_{\lambda}$, due to a deficit of data relative to the background predictions at high values of the BDT score. The $b\bar{b}\gamma\gamma$ limit shows a weaker dependence on $\kappa_{\lambda}$ than the $b\bar{b}b\bar{b}$ and $b\bar{b}\tau^{+}\tau^{-}$ limits because the $b\bar{b}\gamma\gamma$ acceptance varies less as function of $\kappa_{\lambda}$.

The 95% CL allowed $\kappa_{\lambda}$ intervals are given in Table 2. The systematic uncertainties weaken the $\kappa_{\lambda}$ limits by less than 10% relative to those obtained with only statistical uncertainties. The final state least (most) affected by systematic uncertainties is $b\bar{b}\gamma\gamma$ ($b\bar{b}b\bar{b}$). The Higgs boson branching fraction depends on $\kappa_{\lambda}$ due to NLO electroweak corrections [20]. This dependence is neglected in the present treatment, but its overall impact on the allowed $\kappa_{\lambda}$ interval is evaluated to be no more than 7%. Theory uncertainties on the signal cross section shown in Fgure 4(b) are not taken into account when computing the $\kappa_{\lambda}$ limits in Table 2, they affect the limit by less than 8%.

6 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_{\text{KK}}$.

The search for the heavy scalar particle $S$ is performed with all six final states included in this combination. With the exception of $b\bar{b}\tau^{+}\tau^{-}$ and $b\bar{b}b\bar{b}$, 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 $b\bar{b}\tau^{+}\tau^{-}$ 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 $b\bar{b}b\bar{b}$ final state uses the same LO event generator but interfaced to Herwig++ 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$\sigma$ for the search mass range. The combined upper limit on the cross-section is shown as a function of the resonance mass in Figure 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 uncertaintes 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 $b\bar{b}b\bar{b}$, $b\bar{b}W^{+}W^{-}$ and $b\bar{b}\tau^{+}\tau^{-}$ final states only. Gravitons were simulated using an LO model in MadGraph5_aMC@NLO with the NNPDF 2.3 LO PDF set interfaced to Pythia 8.2 with the A14 tune. The resonance width changes with the graviton mass and depends on the parameter $k/\overline{M}_{\text{Pl}}$, where $k$ is the curvature of the warped extra dimension in the bulk RS model and $\overline{M}_{\text{Pl}}=2.4\times 10^{18}\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$ is the effective four-dimensional Planck mass. The search is performed for models with $k/\overline{M}_{\text{Pl}}$ equal to 1 and 2. For $k/\overline{M}_{\text{Pl}}$ = 1 (2), the width ranges from 3% (11%) for a 0.3 TeV graviton mass to 6% (25%) for a 3 TeV graviton mass.

The upper limits in the $G_{\text{KK}}$ search are shown as a function of the resonance mass in Figures 5(b) and 5(c) for $k/\overline{M}_{\text{Pl}}$ equal to 1 and 2, respectively. In the combination, the largest observed deviation from the background expectation is 1.5$\sigma$ (0.7$\sigma$) for the search mass range with $k/\overline{M}_{\text{Pl}}=$1 (2). Exclusion ranges on the $KK$ graviton mass are obtained by comparing the upper limit with the production cross section calculated at LO. In the case of $k/\overline{M}_{\text{Pl}}=1$, the bulk RS model is excluded at 95% CL in the graviton mass range from $306.399$ GeV to $1381.17$ GeV. In the case of $k/\overline{M}_{\text{Pl}}=2$, the model is excluded at 95% CL for graviton masses from 260 GeV, where the scan starts, to $1758.71$ GeV.

The impact of the systematic uncertainties on the upper limits on $G_{\text{KK}}$ has a small dependence on the resonance mass. It is $\sim$20% over the whole mass range for $k/\overline{M}_{\text{Pl}}=1$, and 29% (25%) at a mass of 1 TeV (3 TeV) for $k/\overline{M}_{\text{Pl}}=2$. The largest systematic uncertainties are from the modelling of the backgrounds, affecting the limit by 11% (15%) at 1 TeV (3 TeV) for $k/\overline{M}_{\text{Pl}}=1$ and 16% (21%) at 1 TeV (3 TeV) for $k/\overline{M}_{\text{Pl}}=2$. For $k/\overline{M}_{\text{Pl}}=1$, the subleading systematic uncertainties come from $b$-tagging at low $G_{\text{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 $k/\overline{M}_{\text{Pl}}=2$, subleading systematic uncertainties are from jet energy scale and resolution, impacting the upper limits by 5% at 1 TeV and 4% at 3 TeV. The systematic uncertainties affect upper limits more for $k/\overline{M}_{\text{Pl}}=2$ than for $k/\overline{M}_{\text{Pl}}=1$, because the natural width of the signal graviton is four times larger with $k/\overline{M}_{\text{Pl}}=2$.

7 Constraints on the hMSSM and EWK-singlet models

Exclusion limits are also presented for two specific models, namely the EWK-singlet model [29, 30, 31, 11] and the hMSSM model [26, 69, 27, 28, 11]. The sensitivity of the $b\bar{b}W^{+}W^{-}$, $W^{+}W^{-}W^{+}W^{-}$ and $W^{+}W^{-}\gamma\gamma$ final states to these models is negligible, so the presented results combine only the $b\bar{b}b\bar{b}$, $b\bar{b}\tau^{+}\tau^{-}$ and $b\bar{b}\gamma\gamma$ 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\alpha$ plane (where $m_{S}$ is the resonance mass) for $\tan\beta=1$ and $\tan\beta=2$, shown in Figure 6(a) and Figure 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 ($\Gamma_{S}$) is smaller than the experimental resolution at the resonance mass. This holds when $\Gamma_{S}/m_{S}<2$% for $b\bar{b}\gamma\gamma$, $\Gamma_{S}/m_{S}<5$% for $b\bar{b}b\bar{b}$ and $\Gamma_{S}/m_{S}<10$% for $b\bar{b}\tau^{+}\tau^{-}$. Therefore, the excluded region in the plot is obtained by combining the three final states for $\Gamma_{S}/m_{S}<2$%, by combining the $b\bar{b}b\bar{b}$ and $b\bar{b}\tau^{+}\tau^{-}$ final states for 2% $<\Gamma_{S}/m_{S}<5$%, and using only $b\bar{b}\tau^{+}\tau^{-}$ for 5% $<\Gamma_{S}/m_{S}<10$%. The hatched region shows points where $\Gamma_{S}/m_{S}\geq 10$%, where no exclusion can be provided. Figure 7(a) shows limits for the (sin $\alpha$, tan $\beta$) 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\beta$ plane of the hMSSM model in Figure 7(b). The expected cross-section for each point in the parameter space is obtained using the gluon-gluon fusion cross-section 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\beta$ 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\beta$. The kink at low $\tan\beta$ and high $m_{A}$ values is caused by removing the $b\bar{b}\gamma\gamma$ 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 $b\bar{b}\gamma\gamma$ analysis.

8 Conclusion

A statistical combination of six final states $b\bar{b}b\bar{b}$, $b\bar{b}W^{+}W^{-}$, $b\bar{b}\tau^{+}\tau^{-}$, $W^{+}W^{-}W^{+}W^{-}$, $b\bar{b}\gamma\gamma$ and $W^{+}W^{-}\gamma\gamma$, is presented for the search for non-resonant and resonant production of Higgs boson pairs. These searches use up to 36.1 ${\mathrm{fb}^{-1}}$ of proton–proton collision data at 13 TeV recorded with the ATLAS detector at the LHC111All results are available in digital format on HEPDATA at the following link: https://www.hepdata.net/record/90521.

In both resonant and non-resonant searches, no statistically significant excess of events above the Standard Model predictions is found. For the Standard Model $HH$ production mode, the observed (expected) 95% confidence level upper limit on the gluon–gluon fusion $pp\to HH$ cross-section is $6.91438$ ($10.2769$) times the Standard Model prediction. The expected limit is comparable to the CMS result, while the observed limit is significantly stronger than CMS’s due to a data deficit compared to expected background in ATLAS and an excess in CMS. For the resonant case, upper limits are set on the production cross-section of heavy spin-0 and spin-2 resonances decaying into pairs of Higgs bosons in the mass range 260–3000 GeV.

Upper limits on the $pp\to HH$ cross-section are also computed as a function of the Higgs boson self-coupling modifier $\kappa_{\lambda}=\lambda_{HHH}/\lambda_{HHH}^{\text{SM}}$, by combining the $b\bar{b}b\bar{b}$, $b\bar{b}\tau^{+}\tau^{-}$ and $b\bar{b}\gamma\gamma$ final states. The combination excludes $\kappa_{\lambda}$ values outside the range -4.97183$<\kappa_{\lambda}<$12.0122 (-5.81589$<\kappa_{\lambda}<$12.0083) at 95% confidence level in observation (expectation). The three final states are also combined to constrain the Electroweak Singlet Model in the ($m_{S}$, sin $\alpha$) and the (sin $\alpha$, tan $\beta$) parameter spaces and the habemus Minimal Supersymmetric Standard Model in the ($m_{A}$, tan $\beta$) parameter space.

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, Argentina; YerPhI, Armenia; ARC, Australia; BMWFW and FWF, Austria; ANAS, Azerbaijan; SSTC, Belarus; CNPq and FAPESP, Brazil; NSERC, NRC and CFI, Canada; CERN; CONICYT, Chile; CAS, MOST and NSFC, China; COLCIENCIAS, Colombia; MSMT CR, MPO CR and VSC CR, Czech Republic; DNRF and DNSRC, Denmark; IN2P3-CNRS, CEA-DRF/IRFU, France; SRNSFG, Georgia; BMBF, HGF, and MPG, Germany; GSRT, Greece; RGC, Hong Kong SAR, China; ISF and Benoziyo Center, Israel; INFN, Italy; MEXT and JSPS, Japan; CNRST, Morocco; NWO, Netherlands; RCN, Norway; MNiSW and NCN, Poland; FCT, Portugal; MNE/IFA, Romania; MES of Russia and NRC KI, Russian Federation; JINR; MESTD, Serbia; MSSR, Slovakia; ARRS and MIZŠ, Slovenia; DST/NRF, South Africa; MINECO, Spain; SRC and Wallenberg Foundation, Sweden; SERI, SNSF and Cantons of Bern and Geneva, Switzerland; MOST, Taiwan; TAEK, Turkey; STFC, United Kingdom; DOE and NSF, United States of America. In addition, individual groups and members have received support from BCKDF, CANARIE, CRC and Compute Canada, Canada; COST, ERC, ERDF, Horizon 2020, and Marie Skłodowska-Curie Actions, European Union; Investissements d’ Avenir Labex and Idex, ANR, France; DFG and AvH Foundation, Germany; Herakleitos, Thales and Aristeia programmes co-financed by EU-ESF and the Greek NSRF, Greece; BSF-NSF and GIF, Israel; CERCA Programme Generalitat de Catalunya, Spain; The Royal Society and Leverhulme Trust, United Kingdom.

The crucial computing support from all WLCG partners is acknowledged gratefully, in particular from CERN, the ATLAS Tier-1 facilities at TRIUMF (Canada), NDGF (Denmark, Norway, Sweden), CC-IN2P3 (France), KIT/GridKA (Germany), INFN-CNAF (Italy), NL-T1 (Netherlands), PIC (Spain), ASGC (Taiwan), RAL (UK) and BNL (USA), the Tier-2 facilities worldwide and large non-WLCG resource providers. Major contributors of computing resources are listed in Ref. [75].