Search for non-resonant Higgs boson pair production in the $bb\ell\nu\ell\nu$ final state with the ATLAS detector in $pp$ collisions at $\sqrt{s} = 13$ TeV

TL;DR

This paper reports a search for non-resonant Higgs boson pair production in the $bb\,\ell\nu\ell\nu$ final state using ATLAS data, setting upper limits on the production cross-section relative to the Standard Model prediction.

Contribution

It introduces a neural network-based discriminant to identify Higgs pair production events in the $bb\,\

Findings

Set an upper limit of 1.2 pb on the production cross-section.

Observed data is consistent with the Standard Model within uncertainties.

The analysis improves sensitivity using machine learning techniques.

Abstract

A search for non-resonant Higgs boson pair production, as predicted by the Standard Model, is presented, where one of the Higgs bosons decays via the $H\rightarrow bb$ channel and the other via one of the $H \rightarrow WW^*/ZZ^*/\tau\tau$ channels. The analysis selection requires events to have at least two $b$-tagged jets and exactly two leptons (electrons or muons) with opposite electric charge in the final state. Candidate events consistent with Higgs boson pair production are selected using a multi-class neural network discriminant. The analysis uses 139 fb$^{-1}$ of $pp$ collision data recorded at a centre-of-mass energy of 13 TeV by the ATLAS detector at the Large Hadron Collider. An observed (expected) upper limit of 1.2 ($0.9^{+0.4}_{-0.3}$) pb is set on the non-resonant Higgs boson pair production cross-section at 95% confidence level, which is equivalent to 40 ($29^{+14}_{-9}$) times the value predicted in the Standard Model.

Tables

Table 1. Table 1: List of the ME generators and PS/UE modelling algorithms used in the simulation. Alternative generators and PS/UE models, used to estimate systematic uncertainties, are shown in parentheses. The PDF sets, tunes, and the perturbative QCD highest-order accuracy (leading-order, LO; next-to-leading-order, NLO; next-to-next-to-leading-order, NNLO; next-to-next-to-leading-logarithm, NNLL) used for the normalisation of the samples are also included. The top-quark mass is set to 172.5  GeV .

Process	ME generator	ME PDF	PS/UE model	UE tune	Prediction order for
	(alternative)		(alternative)		total cross-section
$t\overline{t}$ [52, 53]	Powheg-Box v2 [54, 55]	NNPDF3.0NLO [56]	Pythia 8.230 [57]	A14 [58]	NNLO + NNLL [59, 60, 61, 62, 63, 64, 65]
	(MadGraph 5_aMC@NLO)		(Herwig 7.0.4)	(H7-MMHT14)
Single-top $s$-channel, $Wt$ [52, 66, 67]	Powheg-Box	NNPDF3.0NLO	Pythia 8.230	A14	NLO + NNLL [68, 69]
	(MadGraph 5_aMC@NLO)		(Herwig 7.0.4)	(H7-MMHT14)
Single-top $t$-channel [52, 66]	Powheg-Box, MadSpin [48]	NNPDF3.04fNLO	Pythia 8.230	A14	NLO + NNLL [70]
	(MadGraph 5_aMC@NLO)		(Herwig 7.0.4)	(H7-MMHT14)
$W,Z/{\gamma }^{\ast }+\mathrm{jets}$ [71]	Sherpa 2.2.1 [72, 73]	NNPDF3.0NNLO	Sherpa 2.2.1	Sherpa default	NLO(LO) $\le$2(4) partons [74, 75, 76, 77, 78]
($Z/{\gamma }^{\ast }+\mathrm{jets}$)	(MadGraph 5_aMC@NLO)		(Pythia 8.230)	(A14)
Diboson ($WW,WZ,ZZ$) [79]	Sherpa 2.2.2	NNPDF3.0NNLO	Sherpa 2.2.2	Sherpa default	NLO(LO) $\le$1(3) partons [75, 76, 77, 78]
$t\overline{t}W$,$t\overline{t}Z$ [80]	MadGraph 5_aMC@NLO[81]	NNPDF3.0NLO	Pythia 8.210	A14	NLO [82, 83]
$t\overline{t}H$ [80]	MadGraph 5_aMC@NLO	NNPDF3.0NLO	Pythia 8.210	A14	NLO [84, 85]
$WH,ZH$ [86]	Pythia 8.186 [43]	NNPDF2.3LO [45]	Pythia 8.186	A14	NNLO QCD + NLO EW [87, 88, 89, 90, 91, 92, 93]
ggF $H$ [94]	Powheg-Box v2 NNLOPS [95]	CT10 [96]	Pythia 8.212	AZNLO [97]	NNNLO QCD + NLO EW [98]
SM $HH\to bb\ell \nu \ell \nu$ [99]	MadGraph 5_aMC@NLO 2.6.2	CT10	Herwig 7.0.4 [100]	H7-MMHT14 [101]	NNLO [14, 15, 16, 17, 18, 19, 20]

Table 2. Table 2: Description of the variables used as inputs to the DNN classifier.

($p_T$, $\eta$, $\varphi$)	$p_T$, $\eta$, and $\varphi$ of the leptons, leading two jets (not necessarily $b$-tagged), and leading two $b$-tagged jets
Dilepton flavour	Whether the event is composed of two electrons, two muons, or one of each (encoded as 3 booleans)
$\Delta R_{\ell \ell }$, $\left|\Delta {\varphi }_{\ell \ell }\right|$	$\Delta R$ and magnitude of the $\Delta \varphi$ between the two leptons
$m_{\ell \ell }$, $p_T^{\ell \ell }$	Invariant mass and the transverse momentum of the dilepton system
$E_{\text{T}}^{\text{miss}}$, $E_{\text{T}}^{\text{miss}}$-$\varphi$	Magnitude of the missing transverse momentum vector and its $\varphi$ component
$\left|\Delta \varphi ({𝐩}_{\text{T}}^{\text{miss}},{𝐩}_{\text{T}}^{\ell \ell })\right|$	Magnitude of the $\Delta \varphi$ between the ${𝐩}_{\text{T}}^{\text{miss}}$ and the transverse momentum of the dilepton system
$\left|{𝐩}_{\text{T}}^{\text{miss}}+{𝐩}_{\text{T}}^{\ell \ell }\right|$	Magnitude of the vector sum of the ${𝐩}_{\text{T}}^{\text{miss}}$ and the transverse momentum of the dilepton system
Jet multiplicities	Numbers of $b$-tagged and non-$b$-tagged jets
$\left|\Delta {\varphi }_{bb}\right|$	Magnitude of the $\Delta \varphi$ between the leading two $b$-tagged jets
$m_{\text{T}2}^{bb}$	$m_{\text{T}2}$ [120] using the leading two $b$-tagged jets as the visible inputs and ${𝐩}_{\text{T}}^{\text{miss}}$ as invisible input
$H_{\text{T}2}$	Scalar sum of the magnitudes of the momenta of the $H\to \ell \nu \ell \nu$ and $H\to bb$ systems,
	$H_{\text{T}2}=\left|{𝐩}_{\text{T}}^{\text{miss}}+{𝐩}_{\text{T}}^{\ell ,0}+{𝐩}_{\text{T}}^{\ell ,1}\right|+\left|{𝐩}_{\text{T}}^{b,0}+{𝐩}_{\text{T}}^{b,1}\right|$
$H_{\text{T}2}^{\text{R}}$	Ratio of $H_{\text{T}2}$ and scalar sum of the transverse momenta of the $H$ decay products,
	$H_{\text{T}2}^{\text{R}}=H_{\text{T}2}/\left(E_{\text{T}}^{\text{miss}}+\left|{𝐩}_{\text{T}}^{\ell ,0}\right|+\left|{𝐩}_{\text{T}}^{\ell ,1}\right|+\left|{𝐩}_{\text{T}}^{b,0}\right|+\left|{𝐩}_{\text{T}}^{b,1}\right|\right)$,
	where ${𝐩}_{\text{T}}^{\ell \left(b\right),0\left\{1\right\}}$ are the transverse momenta of the leading {subleading} lepton ($b$-tagged jet)

Table 3. Table 3: Analysis region and background estimation summary. Shown are the definitions of the control, validation, and signal regions used in the analysis as well as the predicted and observed event yields in each of these regions. The predicted yields are shown after background-only fits to data in the control regions. The Top and Z / γ ∗ ​ + HF 𝑍 superscript 𝛾 + HF Z/\gamma^{*}\text{+ HF} post-fit normalisation factors, obtained from background-only fits in the corresponding control regions, are shown at the bottom of the table. Also shown is the predicted H ​ H → b ​ b ​ ℓ ​ ν ​ ℓ ​ ν → 𝐻 𝐻 𝑏 𝑏 ℓ 𝜈 ℓ 𝜈 HH\rightarrow bb\ell\nu\ell\nu signal yield in each of the regions, multiplied by a factor of 20. Of the H ​ H 𝐻 𝐻 HH yield in the signal regions, 90% comes from the H ​ H → b ​ b ​ W ​ W ∗ → 𝐻 𝐻 𝑏 𝑏 𝑊 superscript 𝑊 HH\rightarrow bbWW^{*} process, 9% from the H ​ H → b ​ b ​ τ ​ τ → 𝐻 𝐻 𝑏 𝑏 𝜏 𝜏 HH\rightarrow bb\tau\tau process, and 1% from the H ​ H → b ​ b ​ Z ​ Z ∗ → 𝐻 𝐻 𝑏 𝑏 𝑍 superscript 𝑍 HH\rightarrow bbZZ^{*} process. The uncertainties in each yield and in the normalisation correction factors account for the statistical and systematic uncertainties described in Section 6 , with those on the normalisation corrections due only to experimental sources.

Region Definitions
Observable		CR-Top	VR-1	CR-Z+HF	VR-2	SR-SF	SR-DF
Dilepton Flavour		DF	SF	DF or SF	SF	SF	DF
$m_{\ell \ell }$ [GeV]		$(20,60)$	$(20,60)$	$(81.2,101.2)$	$(71.2,81.2)$	$(20,60)$	$(20,60)$
					or $(101.2,115)$
$m_{bb}$ [GeV]		$\notin (100,140)$	$>140$	$(100,140)$	$(100,140)$	$(110,140)$	$(110,140)$
$d_{HH}$		$>4.5$	$>4.5$	$>0$	$>0$	$>5.45$	$>5.55$
Event Yields
Data		108	171	852	157	16	9
Total Bkg.		$108\pm 10$	$162\pm 10$	$852\pm 29$	$147\pm 11$	$14.9\pm 2.1$	$4.9\pm 1.2$
Top		$92\pm 11$	$77\pm 10$	$55\pm 7$	..$71\pm 10$	..$4.8\pm 1.4$	$3.8\pm 1.1$
$Z/{\gamma }^{\ast }\text{+ HF}$		$3.2\pm 0.5$	$70\pm 4$	$686\pm 33$	$60\pm 4$	..$7.8\pm 1.4$	$0.21\pm 0.05$
Other		$13.1\pm 3.4$	$14.2\pm 1.9$	$110\pm 13$	$15.8\pm 1.2$	..$2.3\pm 0.5$	$0.9\pm 0.4$
$HH$ ($\times 20$)		$2.70\pm 0.25$	$1.03\pm 0.22$	.$1.97\pm 0.11$	..$1.22\pm 0.05$	..$5.0\pm 0.6$	$4.8\pm 0.8$
Post-fit Normalisation
${\mu }_{\text{Top}}=0.79\pm 0.10$				${\mu }_{Z/{\gamma }^{\ast }\text{+ HF}}=1.36\pm 0.07$

Table 4. Table 4: Breakdown of the main uncertainty components in the background estimates in the two signal regions for the Top , Z / γ ∗ ​ + HF 𝑍 superscript 𝛾 + HF Z/\gamma^{*}\text{+ HF} , and all other (“Other”) backgrounds. The uncertainty components associated with the total background estimate in the signal regions (the sum of Top , Z / γ ∗ ​ + HF 𝑍 superscript 𝛾 + HF Z/\gamma^{*}\text{+ HF} , and Other) is listed under “Total Bkg.”. As in the upper half of Table  3 , all uncertainties are shown “post-fit”. The percentages show the size of the uncertainty relative to the expected background in each column and uncertainties can be correlated, not necessarily adding in quadrature to the total uncertainty in each column or across each row. Uncertainties in the post-fit normalisation factors, μ Top subscript 𝜇 Top \mu_{\text{Top}} and μ Z / γ ∗ ​ + HF subscript 𝜇 𝑍 superscript 𝛾 + HF \mu_{Z/\gamma^{*}\text{+ HF}} , are only applicable for the Top and Z / γ ∗ ​ + HF 𝑍 superscript 𝛾 + HF Z/\gamma^{*}\text{+ HF} processes.

Uncertainty $[\%]$	SR-SF				SR-DF
	Top	$Z/{\gamma }^{\ast }\text{+ HF}$	Other	Total Bkg.	Top	$Z/{\gamma }^{\ast }\text{+ HF}$	Other	Total Bkg.
Total uncertainty	$28$	$18$	$20$	$14$	$30$	$26$	$41$	$25$
Theoretical	$21$	$15$	$17$	$11$	$20$	$15$	$40$	$17$
Experimental	$12$		$8$		$15$	$17$	$8$	$12$
MC statistics	$8$	$8$	$6$	$8$	$13$	$13$	$7$	$11$
${\mu }_{\text{Top}}$, ${\mu }_{Z/{\gamma }^{\ast }\text{+ HF}}$	$13$	$5$	n/a	$5$	$13$	$5$	n/a	$10$

Table 5. Table 5: Observed and expected upper limits on the ggF-initiated non-resonant H ​ H 𝐻 𝐻 HH production cross-section at 95 % percent 95 95\% CL and their ratios to the SM prediction ( σ SM ​ ( g ​ g → H ​ H ) = 31.05 ± 1.90 superscript 𝜎 SM → 𝑔 𝑔 𝐻 𝐻 plus-or-minus 31.05 1.90 \sigma^{\text{\tiny{SM}}}(gg\rightarrow HH)=31.05\pm 1.90 fb [ 13 , 14 , 15 , 16 , 17 , 18 , 19 , 20 ] ). The ± 1 ​ σ plus-or-minus 1 𝜎 \pm 1\sigma and ± 2 ​ σ plus-or-minus 2 𝜎 \pm 2\sigma variations about the expected limit are also shown. Uncertainties in the SM cross-section are taken into account when computing the upper limits on the cross-section ratio.

	$-2\sigma$	$-1\sigma$	Expected	$+1\sigma$	$+2\sigma$	Observed
$\sigma \left(gg\to HH\right)$ [pb]	0.5	0.6	0.9	1.3	1.9	1.2
$\sigma \left(gg\to HH\right)/{\sigma }^{\text{SM}}\left(gg\to HH\right)$	14	20	29	43	62	40

Full text

\AtlasTitle

Search for non-resonant Higgs boson pair production in the $bb\ell\nu\ell\nu$ final state with the ATLAS detector in $pp$ collisions at $\sqrt{s}=13$ TeV \PreprintIdNumberCERN-EP-2019-143 \AtlasAbstract A search for non-resonant Higgs boson pair production, as predicted by the Standard Model, is presented, where one of the Higgs bosons decays via the $H\rightarrow bb$ channel and the other via one of the $H\rightarrow WW^{*}/ZZ^{*}/\tau\tau$ channels. The analysis selection requires events to have at least two $b$-tagged jets and exactly two leptons (electrons or muons) with opposite electric charge in the final state. Candidate events consistent with Higgs boson pair production are selected using a multi-class neural network discriminant. The analysis uses 139 fb*-1* of $pp$ collision data recorded at a centre-of-mass energy of 13 TeV by the ATLAS detector at the Large Hadron Collider. An observed (expected) upper limit of 1.2 ($0.9^{+0.4}_{-0.3}$) pb is set on the non-resonant Higgs boson pair production cross-section at 95% confidence level, which is equivalent to 40 ($29^{+14}_{-9}$) times the value predicted in the Standard Model.

\AtlasRefCodeHDBS-2018-33 \AtlasJournalPhys. Lett. B 801 (2020) 135145 \AtlasDOI10.1016/j.physletb.2019.135145

1 Introduction

In 2012, the ATLAS and CMS Collaborations reported the observation of a new particle in the search for the Standard Model (SM) Higgs boson ($H$) [1, 2]. So far, measurements of the spin and couplings of the new particle are consistent with those predicted by the Brout–Englert–Higgs (BEH) mechanism of the SM [3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. The SM predicts non-resonant production of Higgs boson pairs ($HH$) in proton–proton ($pp$) collisions, referred to as non-resonant $HH$ production, with the dominant production modes at the LHC proceeding via the gluon–gluon fusion (ggF) process. The ggF process has two leading order contributions: the first corresponds to the so-called ‘triangle diagram’, which includes the Higgs boson self-coupling, and the second is the so-called ‘box diagram’, which includes a heavy-quark loop with two fermion–fermion–Higgs ($ffH$) vertices. These two amplitudes interfere destructively, resulting in a low cross-section of only $31.05\pm 1.90$ fb for the ggF $HH$ production mode, computed at next-to-next-to-leading order (NNLO) and including finite top-quark mass effects [13, 14, 15, 16, 17, 18, 19, 20]. Feynman diagrams illustrating these two contributions are shown in Figure 1. The measurement of non-resonant $HH$ production at the LHC stands as an important test of the BEH mechanism. In many beyond-the-SM (BSM) theories, $HH$ production can be enhanced by modifying the Higgs boson self-coupling, $\lambda_{HHH}$, or the top-quark Yukawa coupling, $y_{t}$, and/or by introducing new contact interactions between two top-quarks or gluons and two Higgs bosons or introducing production mechanisms via intermediate BSM particles [21, 22, 23].

The ATLAS and CMS Collaborations have performed searches for non-resonant $HH$ production in a variety of final states at 13 TeV  [24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. No significant excess of events beyond SM expectations is observed in these searches, with the ATLAS and CMS data-analyses setting observed (expected) limits on non-resonant $HH$ production to be no larger than 6.9 (10.0) and 22.2 (12.8) times the predicted rate in the SM, respectively [34, 35].

This Letter describes a search for non-resonant $HH$ production in the $bb\ell\nu\ell\nu$ final state, where $\ell$ refers to a lepton (either an electron or a muon), using 13 TeV $pp$ collision data collected with the ATLAS detector during 2015–2018 and corresponding to a total integrated luminosity of 139 fb*-1*. The analysis uses machine-learning techniques based on feedforward neural network architectures [36] to construct an event-level classifier trained to distinguish between the $HH$ signal and SM backgrounds. Analyses searching for non-resonant $HH$ production via similar decay channels were performed previously in the single-lepton final state by ATLAS in searches for $HH\rightarrow bbWW^{*}$ [28] and in the dilepton channel by CMS in searches for $HH\rightarrow bbWW^{*}/bbZZ^{*}$ [31].

2 ATLAS detector

The ATLAS detector [37, 38, 39] is a general-purpose particle detector with forward–backward symmetric cylindrical geometry.111 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the $z$-axis along the beam pipe. The $x$-axis points from the IP to the centre of the LHC ring, and the $y$-axis points upwards. Cylindrical coordinates $(r,\phi)$ are used in the transverse plane, $\phi$ being the azimuthal angle around the $z$-axis. The pseudorapidity is defined in terms of the polar angle $\theta$ as $\eta=-\ln\tan(\theta/2)$. The angular distance is measured in units of $\Delta R\equiv\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}$. It includes an inner tracking detector (ID), immersed in an axial magnetic field, which provides precision tracking of charged particles over the range of $|\eta|<2.5$. Calorimeter systems with either liquid argon or scintillator tiles as the active medium provide energy measurements over the range of $|\eta|<4.9$. The muon spectrometer (MS) is positioned outside the calorimeters and includes three air-core toroidal magnets. The MS is composed of several types of muon detectors which provide trigger and high-precision tracking capabilities for $|\eta|<2.4$ and $|\eta|<2.7$, respectively. A hardware-based trigger followed by a software-based trigger reduce the recorded event rate to an average of 1 kHz [40].

3 Dataset and simulated events

The data used for this search were collected in $pp$ collisions at the LHC with a centre-of-mass energy of 13 TeV. Only those data collected during stable LHC beam conditions and with all ATLAS detector subsystems fully operational are used, and correspond to an integrated luminosity of 139 fb*-1*. The selection of candidate events with oppositely charged leptons is based on a combination of single-lepton and dilepton triggers.222Distinct sets of single-lepton triggers are used for electrons and muons. Dilepton triggers require either two electrons, two muons, or one electron and one muon. The use of a given trigger depends on the flavour and the transverse momenta ($p_{\text{T}}$) of the two ($p_{\text{T}}$-ordered) leptons in the event, and on the data-taking period. Single-lepton triggers with $p_{\text{T}}$ thresholds between $22$ and $28\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$ are given priority over dilepton triggers. The criteria of the dilepton triggers are checked only if no single-lepton trigger criteria are met and have $p_{\text{T}}$ thresholds as low as 19 (10)  GeV for the leading (subleading) lepton. At least one reconstructed lepton (or lepton pair) has to match a corresponding trigger object, in which case their offline $p_{\text{T}}$ must be higher than the trigger threshold by at least $2$ GeV, in order to be on the efficiency plateau of the corresponding trigger.

Monte Carlo (MC) simulation [41] is used to model the signal processes and in the estimation of SM background processes. A GEANT4 [42] simulation of the ATLAS detector was used for the background processes. The signal MC samples were processed with a fast simulation that relies on a parameterisation of the calorimeter response [41] and on GEANT4 for the tracking detectors. Simulated events are reconstructed using the same algorithms as used for data and include the effects of multiple $pp$ interactions in the same or neighbouring bunch crossings, collectively referred to as pile-up. The simulation of pile-up collisions was performed with Pythia 8.186 [43] using the ATLAS A3 set of tuned parameters [44] and the NNPDF2.3LO parton distribution function (PDF) set [45]. Simulated events were reweighted to match the distribution of pile-up interactions in data. The average amount of pile-up in the data collected during 2015–2018 was $33.7$.

The signal processes with ggF-initiated non-resonant $HH$ production in the $bb\ell\nu\ell\nu$ final state were generated with an effective Lagrangian in the infinite top-quark mass approximation. The generated signal events were reweighted with form factors that take into account the finite mass of the top-quark [46, 47]. SM background processes were simulated using different MC event generators. The MC matrix element (ME) event generators and PDF sets, the parton showering (PS) and the underlying event (UE) modelling, UE tuned parameters (tune), and the accuracy of the theoretical cross-sections used to normalise the simulated processes are summarised in Table 1. Each SM background process is normalised to the best available respective theoretical cross-section. The mass of the Higgs boson was set to 125 GeV for all signal and background processes. The $HH$ branching fractions (BF) predicted by the SM [13] are used for all Higgs boson decays. MadSpin [48] was used to model top-quark spin correlations and EvtGen [49] was used to model properties of $b$- and $c$-hadron decays for processes using Pythia and for the signal processes.

SM top-quark pair production ($t\bar{t}$) and the production of single top-quarks in association with $W$ bosons ($Wt$) contribute with significant background contamination in the $bb\ell\nu\ell\nu$ final state. At next-to-leading-order (NLO) accuracy, there exists non-trivial interference between these two processes that may be enhanced in phase-space regions wherein there are high fractions of $Wt$ events [50]. Two schemes are typically used to remove the overlap between these two processes: the so-called diagram removal (DR) and diagram subtraction (DS) schemes [51]; the former is used in the present analysis to remove the overlapping events and the latter is used to evaluate the systematic uncertainty in corresponding background event yields. Because of these effects, the sum of the simulated $t\bar{t}$ and $Wt$ processes is considered as a single background process and referred to as the ‘Top’ process in what follows.

4 Event selection and object definitions

Selected events are required to have at least one $pp$ interaction vertex reconstructed from at least two ID tracks with $p_{\text{T}}>0.4\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$. The primary vertex for each event is defined as the vertex with the highest $\sum\left(p_{\text{T}}\right)^{2}$ of associated ID tracks [102]. Events that contain at least one jet arising from non-collision sources or detector noise are rejected by a set of quality criteria [103].

Loose and signal criteria are defined in order to select reconstructed lepton and jet candidates, where the latter is a subset of the former. Compared to the loose objects, the signal objects are required to satisfy tighter identification or quality criteria that are designed to suppress background contributions. Reconstructed loose (signal) electrons are required to satisfy the ‘Loose’ (‘Tight’) likelihood identification criteria [104]. Loose electrons are required to have $p_{\text{T}}>10\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$ and to be within $|\eta|<2.47$. In addition, signal electrons are required to be outside the range $1.37<|\eta|<1.52$, which corresponds to the transition regions between the barrel and endcaps of the electromagnetic calorimeters. In order to reduce background contributions from jets misidentified as electrons, signal electrons are required to be isolated according to the ‘Gradient’ selection criteria [104]. Reconstructed loose and signal muon candidates are required to have $p_{\text{T}}>10\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$, to be within $|\eta|<2.4$, and to satisfy the ‘Medium’ identification criteria [105]. Additionally, signal muons are required to be isolated according to the ‘FixedCutLoose’ selection criteria [105]. Signal electron (muon) candidates are required to originate from the primary vertex by demanding that the significance of the transverse impact parameter, defined as the absolute value of the track transverse impact parameter, $d_{0}$, measured relative to the primary vertex, divided by its uncertainty, $\sigma_{d_{0}}$, satisfy $|d_{0}|/\sigma_{d_{0}}<5$ $(3)$. The difference $\Delta z_{0}$ between the value of the $z$ coordinate of the point on the track at which $d_{0}$ is defined and the longitudinal position of the primary vertex is required to satisfy $|\Delta z_{0}\times\sin\theta|<0.5$ mm, where $\theta$ is the polar angle of the track with respect to the $z$-axis.

Jets are reconstructed from topological clusters of energy deposits in the calorimeters [106] using the anti-$k_{t}$ algorithm [107, 108] with a radius parameter of $R=0.4$ and calibrated as described in Ref. [109]. Candidate loose jets are required to have $p_{\text{T}}>20\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$. Signal jets are required to have $|\eta|<2.8$ and must satisfy pile-up suppression requirements based on the output of a multivariate classifier [110], which identifies jets consistent with a primary vertex in the region $|\eta|<2.4$ and $p_{\text{T}}<120\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$. The MV2c10 multivariate algorithm [111] is used to identify jets containing $b$-hadrons ($b$-tagged jets). An MV2c10 working point with a $b$-tagging efficiency of $70\%$, estimated from simulated $t\bar{t}$ events [112], is used. The $b$-tagged jets must have $p_{\text{T}}>20\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$ and $|\eta|<2.5$. The momentum of $b$-tagged jets is adjusted using the muon-in-jet correction, as described in Ref. [6], by accounting for momentum losses due to muons originating from in-flight semileptonic $b$-hadron decays occurring within the $b$-tagged jet.

The missing transverse momentum $\mathbf{p}_{\text{T}}^{\text{miss}}$, the magnitude of which is denoted by $E_{\text{T}}^{\text{miss}}$, is constructed from the negative vectorial sum of the transverse momenta of calibrated loose objects in the event. An additional term is included to account for the energy of ID tracks that are matched to the primary vertex in the event but not to any of the selected loose objects [113].

To avoid double-counting, loose objects are subject to the overlap removal procedure defined as follows. If a reconstructed electron and muon share a track in the ID, the electron is removed. However, if the muon sharing the track with the electron is calorimeter-tagged,333A calorimeter-tagged muon has only a reconstructed track in the ID matched to energy deposits in the calorimeter compatible with a minimum ionising particle, but no corresponding track segment in the MS. then the muon is removed instead of the electron. If a jet and an electron are reconstructed within $\Delta R=0.2$ of each other, then the jet is removed. If a jet and a muon are within $\Delta R=0.2$ of each other, and the jet has less than three tracks or carries less than $50\%$ of the muon $p_{\text{T}}$, then the jet is removed; otherwise, the muon is removed. Electrons or muons separated from the remaining jets by $\Delta R<0.4$ are removed.

The analysis selects candidate events with exactly two oppositely charged signal leptons, electrons or muons, and at least two signal $b$-tagged jets. To enhance sensitivity to the signal process and to maximise rejection of the expected SM backgrounds, the analysis uses a multivariate approach to select signal events.

5 Analysis strategy

The analysis relies on the use of a multivariate discriminant designed to select candidate events consistent with non-resonant $HH$ production. Section 5.1 describes the architecture and the training of the deep neural network (DNN) classifier from which the discriminant is constructed. Section 5.2 describes the signal region selection criteria. Section 5.3 describes the final background estimation procedure.

5.1 Deep learning approach to target $HH$

The discriminant uses the outputs of a DNN classifier that is built using the Keras library with Tensorflow as a backend [114, 115] and uses the lwtnn library [116] to interface with the analysis software infrastructure of the ATLAS experiment. The sample of events used for training is composed of equal numbers of events from the signal and each of the dominant background processes: Top (as defined in Section 3), $Z/\gamma^{*}\rightarrow\ell\ell$ (Z-$\ell\ell$), and $Z/\gamma^{*}\rightarrow\tau\tau$ (Z-$\tau\tau$) production. The signal sample used in the training of the classifier contains only the $HH\rightarrow bbWW^{*}$ component due to its larger BF relative to the $HH\rightarrow bb\tau\tau$ and $HH\rightarrow bbZZ^{*}$ components. However, the sum of all three signal components is evaluated as the signal when performing the statistical analysis. Additionally, all processes that make up the training sample ($HH\rightarrow bbWW^{*}$, Top, Z-$\ell\ell$, and Z-$\tau\tau$) have the same weight during the training of the classifier. The training sample is composed of simulated candidate events with $m_{\ell\ell}>20\leavevmode\nobreak\ \text{Ge\kern-1.00006ptV}$ and having one or more $b$-tagged jets, where events with exactly one $b$-tagged jet are included to increase the number of events available for training. For the training events with exactly one $b$-tagged jet, each observable that requires at least two $b$-tagged jets is set to its mean value as computed with the full set of training events that contain at least two $b$-tagged jets. Observables that require two $b$-tagged jets are defined using the leading two $b$-tagged jets. The classifier contains two fully connected hidden layers each with 250 nodes. Rectified linear unit (ReLU) activations are used for each layer [117]. In order to improve the robustness of the training and to reduce effects due to overtraining, there is a dropout layer that randomly drops 50% of the nodes between the two fully connected layers during training [118]. The classifier produces four outputs that are passed through a softmax activation, constraining their sum to one [36]. The resulting four outputs, each constrained to values between 0 and 1, are referred to as $p_{i}$ ($i\in\{HH,\text{Top},\text{Z-}\ell\ell,\text{Z-}\tau\tau\}$). Values of $p_{i}$ nearer to 1 indicate that the event likely belongs to class $i$ and values nearer to 0 indicate otherwise. The main discriminant in the analysis, $d_{HH}$, is constructed from the four $p_{i}$ and is defined as $d_{HH}=\ln\left[p_{HH}/\left(p_{\text{Top}}+p_{\text{Z-}\ell\ell}+p_{\text{Z-}\tau\tau}\right)\right]$.

The $HH\rightarrow bb\ell\nu\ell\nu$ signal events are characterised by two distinct ‘Higgs hemispheres’. One hemisphere contains the two $b$-tagged jets from the $H\rightarrow bb$ decay and it is typically opposite in the transverse plane to the second hemisphere that contains the two leptons and $E_{\text{T}}^{\text{miss}}$ from the $H\rightarrow WW^{*}/ZZ^{*}/\tau\tau$ decay. The final-state objects in the SM backgrounds, the Top process in particular, are distributed more uniformly within the event and they typically do not exhibit the same opposite hemispheres topology as the $HH$ signal. These Higgs hemispheres thus provide a topological criterion that distinguishes the signal from the background and motivates the choice of input observables that are provided to the classifier. Thirty-five such variables are provided as inputs to the classifier, ranging from momentum components of the visible final-state objects to observables using event-wide information, and are constructed using only calibrated final state objects that have well-understood uncertainties (Section 6). A complete list is provided in Table 2. The event-wide input observables are sensitive to the presence of Higgs hemispheres in the signal and are largely angular in nature or take advantage of the fact that the final state objects from each of the Higgs bosons in the signal tend to be near to each other. The observables $H_{\text{T}2}^{\text{R}}$ and $m_{\text{T}2}^{bb}$ are non-standard high-level observables that are not straightforward functions of the momenta of the final-state objects. By construction, $H_{\text{T}2}^{\text{R}}$ can take values between zero and one; it peaks near one for signal and is more broadly distributed for background. The $m_{\text{T}2}^{bb}$ observable is defined similarly to the $M_{\text{T}2}^{bb}$ observable in Ref. [119] but does not include the final-state leptons. As discussed in Ref. [119], for the Top backgrounds $m_{\text{T}2}^{bb}$ generally has values below the mass of the top-quark due to kinematic constraints while for the $Z/\gamma^{*}$ processes, which have little-to-no $E_{\text{T}}^{\text{miss}}$, $m_{\text{T}2}^{bb}$ is typically below $45$ GeV. The use of dropout regularisation during the training of the classifier allows it to more effectively use the information contained in the full set of inputs presented in Table 2 by reducing its susceptibility to overtraining effects that may otherwise appear as a result of using such an extended input feature space in the case where no such regularisation is performed. To verify this, the performance of the classifier was checked using an independent sample of events not used in the training of the classifier and was found to be compatible to its performance when presented with those of the training sample.

5.2 Signal selection criteria

To define signal selection criteria, the analysis relies on the invariant mass of the two leptons, $m_{\ell\ell}$, and the invariant mass of the two leading ($p_{\text{T}}$-ordered) $b$-tagged jets, $m_{bb}$. Due to spin-correlation effects present in the $H\rightarrow WW^{*}\rightarrow\ell\nu\ell\nu$ decay within the dominant $HH\rightarrow bbWW^{*}$ signal process, the signal events exhibit values of $m_{\ell\ell}$ that are typically below $60$ GeV. By selecting low values of $m_{\ell\ell}$, the signal purity can therefore be enhanced while rejecting a large component of the SM $Z$ boson and Top backgrounds. Additionally, $m_{bb}$ has a peak at the mass of the Higgs boson for the signal process and therefore provides an effective means to define selections in which the $HH$ contribution is enhanced. The signal selection criteria therefore require $m_{\ell\ell}\in(20,60)$ GeV and $m_{bb}\in(110,140)$ GeV. The $m_{\ell\ell}>20$ GeV requirement is enforced in order to remove contamination from low-mass resonances and $Z/\gamma^{*}$ processes. The signal selection criteria are further broken down into same-flavour (SF), i.e. $ee$ or $\mu\mu$, or different-flavour (DF), i.e. $e\mu$, regions. Separating by dilepton flavour enhances the separation power between the signal and $Z/\gamma^{*}$ background; the former has roughly equal probabilities for the SF and DF final states and the latter leads predominantly to SF final states.

In addition to the $m_{\ell\ell}$ and $m_{bb}$ requirements, the same-flavour and different-flavour signal regions, SR-SF and SR-DF, respectively, are defined by requiring high values of $d_{HH}$ and are presented in Table 3. The chosen threshold values of $d_{HH}>5.45$ ($5.55$) for SR-SF (SR-DF) are found to maximise the expected sensitivity to the non-resonant $HH$ process. The predicted $HH\rightarrow bb\ell\nu\ell\nu$ signal yields in SR-SF and SR-DF are shown in Table 3, and are composed of $90\%$ $HH\rightarrow bbWW^{*}$, $9\%$ $HH\rightarrow bb\tau\tau$, and $1\%$ $HH\rightarrow bbZZ^{*}$. The predominance of the $HH\rightarrow bbWW^{*}$ process over the other two is a result of both its larger overall BF and of the classifier having been trained only on this component of the signal.

5.3 Background estimation

As mentioned in Section 5.1, the dominant backgrounds expected to contaminate the signal regions are the Top and $Z/\gamma^{*}$ processes, specifically $Z/\gamma^{*}$ production in association with jets originating from heavy-flavour hadrons ($bb$, $bc$, or $cc$), subsequently referred to as $Z/\gamma^{*}\text{+ HF}$. Subdominant SM processes contribute via $t\bar{t}$ production in association with an electroweak vector boson, single Higgs boson production (predominantly via the $t\bar{t}H$ mode), $Z/\gamma^{*}$ production in association with light-flavour jets, and electroweak diboson processes. There is additionally a minor contribution of background events from non-prompt leptons produced in semileptonic decays of heavy-flavour hadrons and from misidentified electron candidates arising from photon conversions and jets. This background is estimated using events with a same-charge lepton pair, following procedures described in Ref. [121], after subtracting the prompt-lepton contribution. The rest of the SM background processes detailed in Table 1 are estimated primarily using simulation.

Dedicated control regions are defined to derive data-driven normalisation corrections for the dominant background processes: CR-Top for Top and CR-Z+HF for $Z/\gamma^{*}\text{+ HF}$. These normalisation corrections have a uniform prior and are checked in two validation regions, VR-1 and VR-2, enriched with events from the Top and $Z/\gamma^{*}\text{+ HF}$ processes. The control and validation regions are defined in Table 3 and are kinematically close to the signal regions. CR-Top (CR-Z+HF) and VR-1 (VR-2) are defined by inverting the $m_{bb}$ ($m_{\ell\ell}$) requirements relative to those of the signal regions but retain a selection of the high $d_{HH}$ region similar to the signal regions. The $d_{HH}$ selections were relaxed to increase statistical power, independent checks showed that this did not have a significant impact on the post-fit normalisation corrections in Table 3.

VR-1 keeps only those events with $m_{bb}>140$ GeV, excluding the region $m_{bb}<100$ GeV which is included in CR-Top, due to significant contamination of $Z/\gamma^{*}\text{+ HF}$ events. The correlations between the $m_{\ell\ell}$ and $m_{bb}$ observables and $d_{HH}$ after the preselection are observed to be small and do not prevent the use of the former two in the construction of the analysis regions defined in Table 3, as $d_{HH}$ is found to rely mainly on the information provided by the additional input observables listed in Table 2. This absence of strong correlation ensures that the measurements made in the tails of $d_{HH}$ in the control regions can be extrapolated to those in the signal regions.

The Top background in the signal regions is expected to be composed of approximately equal contributions from the $t\bar{t}$ and single-top-quark $Wt$ process and therefore susceptible to the interference effects as described in Section 3. For this reason, CR-Top and the validation regions are defined so that they have predicted $t\bar{t}$ and $Wt$ compositions similar to that of the signal regions. This ensures that the normalisation correction determined in the fit for the Top background results in an accurate estimate of the combined $t\bar{t}$ and $Wt$ process in the signal regions, accounting for potential interference effects present in data but not necessarily modelled in MC simulation. Table 3 compares the observed and predicted event yields, where the background event yields obtained after background-only fits in the corresponding control regions are also shown. The post-fit normalisation correction factors for the Top and $Z/\gamma^{*}\text{+ HF}$ background processes, respectively $\mu_{\text{Top}}=0.79\pm 0.10$ and $\mu_{Z/\gamma^{*}\text{+ HF}}=1.36\pm 0.07$, are also shown in Table 3. The uncertainties in $\mu_{\text{Top}}$ and $\mu_{Z/\gamma^{*}\text{+ HF}}$ take into account the statistical and systematic uncertainties due to the experimental sources, as described in Section 6.

Distributions of $d_{HH}$ in the control regions after performing background-only fits to data in the control regions and applying the Top and $Z/\gamma^{*}\text{+ HF}$ normalisation corrections are shown in Figure 2. In the control and validation regions, good agreement between the data and SM prediction provided by the post-fit MC simulation is observed for the observables relevant to the analysis.

6 Systematic uncertainties

The analysis evaluates several sources of systematic uncertainty for the signal and background processes, which are classified as either experimental (detector or luminosity related) or theoretical modelling uncertainties. Statistical uncertainties of the simulated event samples are also taken into account. The main uncertainty components are summarised in Table 4. MC modelling uncertainties in the Top and $Z/\gamma^{*}\text{+ HF}$ background estimates are dominant, followed by statistical and detector uncertainties.

The normalisation corrections of the Top and $Z/\gamma^{*}\text{+ HF}$ background processes are determined primarily by the data events in the control regions when performing the statistical analysis. These corrections take into account the statistical and systematic uncertainties due to the experimental sources, as described later in this section. In addition, the systematic uncertainties in the theoretical modelling of these processes are applied as uncertainties in the corrected predictions in the signal regions using the following procedures. The uncertainties in the estimated Top background event yields due to parton shower modelling are assessed as the difference between the predictions of Powheg-Box showered with Pythia or Herwig, and those due to the choice of event generator are assessed by comparing the predictions of Powheg-Box or MadGraph 5_aMC@NLO [122], both showered with Pythia. The uncertainties due to missing higher-order corrections are estimated by changing the renormalisation and factorisation scales ($\mu_{\textsc{r}}$ and $\mu_{\textsc{f}}$, respectively) up and down by a factor of two (8-points variation). The uncertainties due to the modelling of initial- and final-state radiation (ISR and FSR, respectively) in the generators used to simulate the Top background processes are evaluated using the method described in Ref. [122]. The Top background composition is varied within the uncertainties in the theoretical predictions for the $t\bar{t}$ and single-top-quark $Wt$ cross-sections [65, 123, 68]. The uncertainty arising from the interference between the NLO predictions for $t\bar{t}$ and $Wt$ processes is estimated by taking the difference between the predicted Top background yields obtained with the DR and DS schemes used for the NLO $Wt$ calculation [122]. The uncertainties due to PDF variations are computed as the envelope of the central values of the nominal NNPDF3.0 PDF set and the CT14, MMHT14, and PDF4LHC15_30PDF PDF sets [124]. All uncertainties except those in the scale variations, cross-section, and interference are considered as fully correlated between the $t\bar{t}$ and $Wt$ processes. The $Z/\gamma^{*}\text{+ HF}$ modelling uncertainties are estimated using the nominal Sherpa 2.2.1 samples by considering different merging (CKKW-L) [125] and resummation scales. The uncertainties due to PDF variations and changes in $\mu_{\textsc{r}}$ and $\mu_{\textsc{f}}$ are calculated using the same procedures as for the Top backgrounds. An additional uncertainty in the $Z/\gamma^{*}$ process is computed by taking the difference between the nominal Sherpa 2.2.1 samples with samples generated using MadGraph 5_aMC@NLO+Pythia8. The dominant uncertainties in the total background estimates in SR-SF are the $Z/\gamma^{*}\text{+ HF}$ modelling uncertainties ($8\%$), primarily that due to comparison of Sherpa 2.2.1 and MadGraph 5_aMC@NLO, and the parton shower uncertainty affecting the Top background process ($5\%$). The uncertainties in the background estimates in SR-DF are dominated by the uncertainty due to the parton shower affecting the Top background process ($12\%$), the uncertainty in the Top normalisation correction $\mu_{\text{Top}}$ ($10\%$), the uncertainty due to the comparison between the generators used for the Top process ($7.5\%$), and the uncertainty due to the modelling of ISR and FSR in the Top process ($5\%$).

Systematic uncertainties in the signal acceptance due to varying $\mu_{\textsc{r}}$ and $\mu_{\textsc{f}}$, as well as PDF-induced uncertainties, are evaluated using the same procedure as for the Top background process. The resulting scale (PDF) uncertainties are $<3\%$ ($<1\%$) in both signal regions. The uncertainty due to the parton shower modelling is computed by comparing Herwig7 with Pythia8, and is found to be $8\%$ ($9\%$) in SR-SF (SR-DF). The uncertainty in the $HH$ production cross-section, evaluated to be $5\%$, is included as an uncertainty in $\sigma^{\text{\tiny{SM}}}(gg\rightarrow HH)$ when computing the upper limits on the cross-section ratio in Table 5. This value is the quadrature sum of the scale, PDF+$\alpha_{s}$, and top mass contributions as reported by the LHCXSWG [20].

The uncertainties due to experimental sources arise primarily from the mismeasurement of reconstructed object momenta and from the mismodelling of reconstruction efficiencies. These uncertainties include uncertainties from the mismodelling of the jet energy scale (JES) [109] and jet energy resolution (JER) [126]. Additional uncertainties for $b$-tagged jets arise from the mismodelling of the $b$-tagging efficiency [111] and from the mismodelling of the rates at which charm- and light-flavoured jets are selected as $b$-tagged jets [127, 128]. Lepton-related uncertainties arise from the mismodelling of the electron [104] (muon [105]) reconstructed energy (momentum) measurements, as well as in the mismodelling of their reconstruction and identification efficiencies [104, 105]. The $E_{\text{T}}^{\text{miss}}$ scale and resolution [113] uncertainties, as well as uncertainties from the mismodelling of pile-up, trigger efficiency and luminosity, are also taken into account. The uncertainty in the combined 2015–2018 integrated luminosity is $1.7\%$ [129], obtained using the LUCID-2 detector [130] for the primary luminosity measurements. The combined effect of the experimental sources of systematic uncertainty in the predicted background yields is summarised in Table 4 and is dominated by the JER, with all other contributions found to be negligible.

7 Results

In order to extract information about the $HH\rightarrow bb\ell\nu\ell\nu$ signal cross-section, a counting experiment is performed with a profile-likelihood fit [131] simultaneously across the CR-Top, CR-Z+HF, SR-SF, and SR-DF regions using the predicted and observed event counts in each region as inputs. The Top and $Z/\gamma^{*}\text{+ HF}$ normalisation corrections are also extracted from this fit and are found to differ negligibly from those presented in Table 3. All sources of systematic and statistical uncertainty in the signal and background models are implemented as deviations from the nominal model, scaled by nuisance parameters that are profiled in the fit. The $p$-value corresponding to the background-only hypothesis, giving the probability that the data in the signal regions be at least as incompatible with the background-only hypothesis as that observed in SR-SF and SR-DF, is $p_{0}=0.15$ and corresponds to $1.05\sigma$ significance. Distributions of $m_{bb}$, $m_{\ell\ell}$, and $d_{HH}$ after performing background-only fits to data in the control regions and applying the Top and $Z/\gamma^{*}\text{+ HF}$ normalisation corrections are shown in Figure 3. The signal selection criteria are imposed on all observables shown in Figure 3 apart from the one being plotted, except that the $d_{HH}$ requirement for the $m_{bb}$ and $m_{\ell\ell}$ distributions is relaxed to $d_{HH}>5$. No significant excess of events over the expected SM background is observed and upper limits are set on non-resonant Higgs boson pair production at $95\%$ confidence level (CL) using the $\text{CL}_{\text{s}}$ method [132]. Table 5 presents these upper limits and comparisons with the SM prediction. The observed (expected) limit at $95\%$ CL is $1.2$ ($0.9$) pb, corresponding to $40$ ($29$) times the SM prediction.

8 Conclusions

A search for non-resonant Higgs boson pair production, as predicted by the SM, is presented in the final state with at least two $b$-tagged jets and exactly two leptons with opposite electric charge, where one of the Higgs bosons decays to $bb$ and the other decays to either $WW^{*}$, $ZZ^{*}$, or $\tau\tau$. The analysis uses $pp$ collision data recorded at $\sqrt{s}=13$ TeV by the ATLAS detector at the LHC, corresponding to an integrated luminosity of 139 fb*-1*. The data are in agreement with the predictions for the SM background processes. An observed (expected) $95\%$ CL upper limit is set on the cross-section for the production of Higgs boson pairs, corresponding to $40$ ($29$) times the SM prediction. These limits are comparable to the previous leading searches for non-resonant Higgs boson pair production performed by the ATLAS and CMS experiments.

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. [133].