# Search for Higgs boson pair production in the two bottom quarks plus two photons final state in collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

The Atlas Collaboration

## ► To cite this version:

The Atlas Collaboration. Search for Higgs boson pair production in the two bottom quarks plus two photons final state in collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector. 2021. hal-03391686

HAL Id: hal-03391686

https://hal.science/hal-03391686

Preprint submitted on 9 Nov 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Search for Higgs boson pair production in the two bottom quarks plus two photons final state in $pp$ collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

The ATLAS Collaboration

Searches are performed for non-resonant and resonant di-Higgs boson production in the $b\bar{b}\gamma\gamma$ final state. The data set used corresponds to an integrated luminosity of 139 fb$^{-1}$ of proton-proton collisions at a center-of-mass energy of 13 TeV recorded by the ATLAS detector at the CERN Large Hadron Collider. No excess with respect to background expectations is found and upper limits on the di-Higgs boson production cross sections are set. A 95% confidence level upper limit of 130 fb is set on the $pp \to HH$ non-resonant production, where the expected limit is 180 fb. The observed (expected) limit corresponds to 4.1 (5.5) times the cross section predicted by the Standard Model. The observed (expected) limit on the Higgs boson trilinear coupling modifier $\kappa_\lambda$ is extracted to be [-1.5, 6.7] ([-2.4, 7.7]) at 95% confidence level. The constraints on $\kappa_\lambda$ are obtained over an expected hypothesis excluding $pp \to HH$ production. For the resonant production of a new hypothetical scalar particle $X$ ($X \to HH \to b\bar{b}\gamma\gamma$), limits on the cross section $pp \to X \to HH$ are presented for the narrow-width approximation as a function of $m_X$ in the range $251 \text{ GeV} \le m_X \le 1000 \text{ GeV}$. The observed (expected) limits on the cross section $pp \to X \to HH$ range from 610 fb to 47 fb (360–43 fb) over the considered mass range.
# 1 Introduction

Since the discovery of the Higgs boson in 2012 [1, 2], one of the priorities of the ATLAS and CMS Collaborations has been to better understand the properties of the Brout-Englert-Higgs mechanism [3-8]. The Higgs boson self-coupling provides information about the structure of the Higgs potential. A direct probe of the self-coupling of the Higgs boson is possible by studying Higgs boson pair (HH) production. Furthermore, an enhancement of the Higgs boson pair production rate with respect to the Standard Model (SM) prediction would point to new beyond-the-Standard-Model (BSM) physics and may be within the sensitivity reach of the proton-proton (pp) collision data collected at $\sqrt{s} = 13$ TeV during the Run 2 of the Large Hadron Collider (LHC) [9].

At leading order (LO), the production of Higgs boson pairs via gluon-gluon fusion (ggF) proceeds through the two diagrams shown in Figure 1. These diagrams interfere destructively, leading to a small production cross section [10–12]. For 13 TeV pp collisions and a Higgs boson mass $m_H = 125.09$ GeV, the ggF cross section, calculated at next-to-next-to-leading-order (NNLO) accuracy in finite top mass approximation (FTapprox), is $\sigma_{HH}$ (ggF) = $31.02_{-5.0\%}^{+2.2\%}$ (Scale) $\pm$3.0% ($\alpha_s$+PDF) $\pm$2.6% ($m_{top}$) fb [13–15], where ‘Scale’ represents the uncertainty due to missing higher order quantum chromodynamics (QCD) calculations, ‘$\alpha_s$+PDF’ the uncertainties on the strong coupling constant and parton distribution functions, and ‘$m_{top}$’ the uncertainty related to the top-quark mass scheme [16, 17].

<figure><img src="image_2.png" /><figcaption>Figure 1: The Feynman diagrams for the dominant gluon-gluon fusion production processes. In the Standard Model, the (a) trilinear coupling process, (b) the ‘box’ diagram, and the destructive interference between the two processes contribute to the total cross section. In the figure, $\kappa_\lambda$ represents the Higgs boson trilinear coupling modifier, $\kappa_\lambda = \lambda_{HHH} / \lambda_{HHH}^{SM}$.</figcaption></figure>

The di-Higgs vector-boson fusion (VBF) production cross section, calculated at next-to-next-to-next-to-leading-order (N3LO), is $\sigma_{HH}$ (VBF) = $1.72_{-0.04\%}^{+0.03\%}$ (Scale) $\pm$2.1% ($\alpha_s$+PDF) fb [13], for $m_H = 125.09$ GeV, which is one order of magnitude lower than the one from the ggF process. The VBF production mode provides additional sensitivity on the Higgs trilinear coupling, as shown in Figure 2. Both the ggF and VBF production modes of Higgs boson pairs are considered as a signal in this paper. The rest of production modes have lower cross sections [18] and are neglected.

Non-resonant enhancements to the Higgs boson pair cross section can either originate from loop corrections involving new particles, such as light, colored scalars [19], or through non-SM couplings between the Higgs boson and other SM particles. The non-resonant production cross section can also be altered due to the trilinear self-coupling, $\lambda_{HHH}$, being different than the SM prediction, as discussed in Refs. [20, 21]. Such an effect can be captured by a scale factor $\kappa_\lambda$ defined as $\kappa_\lambda = \lambda_{HHH} / \lambda_{HHH}^{SM}$, where $\lambda_{HHH}^{SM}$ is the SM value of the parameter.

In addition to the non-resonant enhancements, searching for resonant production of Higgs boson pairs is well motivated. Figure 3 shows a ggF production diagram possible in BSM theories predicting the
<figure><img src="image_3.png" /><figcaption>Figure 2: The VBF production of Higgs boson pairs via a (a) *HHVV* vertex, (b) trilinear coupling, and (c) *VVH* production mode. In the figure, $\kappa_\lambda$ represents the Higgs boson trilinear coupling modifier, $\kappa_\lambda = \lambda_{HHH}/\lambda_{HHH}^{SM}$.</figcaption></figure>

existence of heavy scalar particles that can decay into a pair of Higgs bosons. Such theories include models with two Higgs doublets [22], like the minimal supersymmetric extension of the SM [23], twin Higgs models [24] and composite Higgs models [25], adding a second complex scalar doublet to the Higgs sector. Alternatively, the Randall–Sundrum model of warped extra dimensions [26] predicts spin-0 radions that could couple to a Higgs boson pair.

<figure><img src="image_4.png" /><figcaption>Figure 3: Gluon-gluon fusion production of a heavy resonance decaying into a Higgs boson pair.</figcaption></figure>

This paper presents a search for di-Higgs production in the $b\bar{b}\gamma\gamma$ final state, including dedicated assessments of non-resonant and resonant contributions. The analysis considers the full Run 2 data set of 139 fb$^{-1}$ at 13 TeV. For both the non-resonant and resonant searches, the analysis employs a multivariate method designed to reject background processes and the statistical results are obtained from a fit of the di-photon invariant mass, $m_{\gamma\gamma}$. For the non-resonant search, data are divided into different categories based on the four-body invariant mass to target different $\kappa_\lambda$ ranges. The resonant search focuses on probing the existence of a narrow-width scalar particle X in the mass range $251 \text{ GeV} < m_X < 1000 \text{ GeV}$ and selection criteria depend on the mass of the probed scalar particle. The main background processes are the di-photon plus jets production and processes where a Higgs boson decaying into a pair of photons is produced. In the context of the resonant search, the non-resonant *HH* production is considered as a background.

Previous results by the ATLAS Collaboration were obtained in this channel with an integrated luminosity of 36 fb$^{-1}$ of data at 13 TeV collected during Run 2, and were found to be consistent with SM expectations within uncertainties [27]. The search for non-resonant enhancements of the Higgs boson pair production set an observed (expected) 95% confidence level (CL) upper limit on the *HH* cross section of 0.73 (0.93) pb, corresponding to 22 (28) times the SM prediction. The Higgs trilinear coupling was constrained to be between -8.2 < $\kappa_\lambda$ < 13.2 at 95% CL (-8.3 < $\kappa_\lambda$ < 13.2 expected). A previous combination of searches for *HH* pair production performed by the ATLAS Collaboration with up to 36 fb$^{-1}$ of 13 TeV data provided constraints of -5.0 < $\kappa_\lambda$ < 12.0 at 95% CL [28]. The search for enhancements due to the decay of a narrow-width scalar particle was performed, and results given as a function of the resonance mass, $m_X$. The
observed (expected) limits were between 1.1 pb and 0.12 pb over the range $260 \text{ GeV} < m_X < 1000 \text{ GeV}$. A combination of searches by the ATLAS Collaboration with up to $36 \text{ fb}^{-1}$ of 13 TeV data for a narrow-width scalar resonance decaying into a $HH$ pair was performed and provided upper limits between 851 fb and 4.6 fb over the range $260 \text{ GeV} < m_X < 3000 \text{ GeV}$ [28]. The CMS Collaboration also set observed (expected) upper limits of 7.7 (5.2) times the SM prediction at 95% CL non-resonant enhancements of the Higgs boson pair production in the $b\bar{b}\gamma\gamma$ final state with $137 \text{ fb}^{-1}$ of 13 TeV data [29]. The CMS Collaboration set 95% CL upper limits on the production cross section of a narrow-width scalar particle between 4.2 pb and 0.23 pb over the range $260 \text{ GeV} < m_X < 900 \text{ GeV}$ with $35.9 \text{ fb}^{-1}$ of 13 TeV data [30].

## 2 The ATLAS detector

The ATLAS detector [31] at the LHC covers nearly the entire solid angle around the collision point.<fnref n="1" /> It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadron calorimeters, and a muon spectrometer incorporating three large superconducting air-core toroidal magnets.

The inner-detector system (ID) is immersed in a 2 T axial magnetic field and provides charged-particle tracking in the range $|\eta| < 2.5$. The high-granularity silicon pixel detector covers the vertex region and typically provides four measurements per track, the first hit normally being in the insertable B-layer (IBL) installed before Run 2 [32, 33]. It is followed by the silicon microstrip tracker (SCT), which usually provides eight measurements per track. These silicon detectors are complemented by the transition radiation tracker (TRT), which enables radially extended track reconstruction up to $|\eta| = 2.0$. The TRT also provides electron identification information.

The calorimeter system covers the pseudorapidity range $|\eta| < 4.9$. Within the region $|\eta| < 3.2$, electromagnetic calorimetry is provided by barrel and endcap high-granularity lead/liquid-argon (LAr) calorimeters, with an additional thin LAr presampler covering $|\eta| < 1.8$ to correct for energy loss in material upstream of the calorimeters. Hadron calorimetry is provided by the steel/scintillator-tile calorimeter, segmented into three barrel structures within $|\eta| < 1.7$, and two copper/LAr hadron endcap calorimeters. The solid angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules optimized for electromagnetic and hadronic energy measurements respectively.

The muon spectrometer (MS) comprises separate trigger and high-precision tracking chambers measuring the deflection of muons in a magnetic field generated by the superconducting air-core toroidal magnets. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. A set of precision chambers covers the region $|\eta| < 2.7$ with three layers of monitored drift tubes, complemented by cathode-strip chambers in the forward region, where the background is highest. The muon trigger system covers the range $|\eta| < 2.4$ with resistive-plate chambers in the barrel, and thin-gap chambers in the endcap regions.

Interesting events are selected by the first-level trigger system implemented in custom hardware, followed by selections made by algorithms implemented in software in the high-level trigger [34]. The first-level

<footnote n="1">¹ ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector and the z-axis along the beam pipe. The x-axis points from the IP to the center 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 rapidity $y$ is defined in terms of the energy, the momentum and the polar angle $\theta$: $y = \frac{1}{2} \ln \left( \frac{E+p \cdot \cos \theta}{E-p \cdot \cos \theta} \right)$. The angular distance is measured in units of $\Delta R \equiv \sqrt{(\Delta\eta)^2 + (\Delta\phi)^2}$.</footnote>
trigger accepts events from the 40 MHz bunch crossings at a rate below 100 kHz, which the high-level
trigger further reduces in order to record events at about 1 kHz rate.

3 Data and simulation samples

This analysis uses pp collision data collected by the ATLAS experiment from 2015 to 2018 with proton
beams colliding at a center-of-mass energy of $\sqrt{s} = 13$ TeV. After data quality requirements [35] the full
data set represents an integrated luminosity of 139.0 $\pm$ 2.4 fb$^{-1}$ [36, 37]. A mean number of 34.2 pp
inelastic interactions per bunch crossing [38] is recorded.

Monte Carlo (MC) simulations are available for the signal as well as most background processes as detailed
in the rest of the section. The reducible backgrounds from final states with jets wrongly identified as
photons (γ-jet and di-jet backgrounds) are however estimated using a data-driven technique detailed in
Section 4.3.

Events from ggF non-resonant *HH* production are generated at next-to-leading-order accuracy in QCD with finite top-quark mass in both real and virtual corrections (NLO FT) [11], using the Powheg Box v2 [39] generator in the finite top-quark mass approximation [40, 41] with the PDF4LHC15 parton distribution function (PDF) set [42]. The Pythia 8.244 generator is used for parton showering, hadronization and underlying event simulation. Herwig v7.1.6 is used as alternative generator to calculate the theory uncertainty from parton shower. Samples are generated for the coupling modifier $\kappa_\lambda = 1$ and 10.

For the ggF non-resonant *HH* production, a reweighting method is used to provide predictions at different
$\kappa_\lambda$ values. The reweighting method derives the scale factors as a function of $\kappa_\lambda$ in bins of $m_{HH}$ by
performing a linear combination of samples generated at different $\kappa_\lambda$ values. This method has been
validated by comparing the $m_{\gamma\gamma}$ distributions and event yields in the generated $\kappa_\lambda = 10$ sample and
the $\kappa_\lambda = 1$ sample reweighted to $\kappa_\lambda = 10$. Good agreement is obtained in all categories. A systematic
uncertainty in the range of 3–4% is associated to the reweighting process, based on the observed maximum
differences. For each $\kappa_\lambda$ value, the inclusive cross section is normalized according to Ref. [43].

Events from the VBF non-resonant *HH* production are generated at LO [44, 45] using
MadGraph5_aMC@NLO v2.6.0 [44]. The NNPDF3.0nlo PDF set [46] is used in the matrix element,
interfaced to Pythia 8.244. The cross section of the VBF *HH* process is evaluated at N3LO in QCD [47–49],
as outlined in Section 1. Since samples are generated at LO for four values of the coupling modifier
$\kappa_\lambda = 0, 1, 2$ and 10, the N3LO to LO cross section ratio at the SM value is calculated and this factor is
applied to the VBF *HH* cross section. These samples are used to derive a parametrization of the signal
yields in the signal region as a function of $\kappa_\lambda$ by performing a fit with a second order polynomial to the MC
predictions in each analysis category to be detailed in Section 4.2.2.

The process of a heavy spin-0 resonance *X* being produced via ggF and decaying into a pair of Higgs
bosons, $pp \to X \to HH$, is simulated using MadGraph5_aMC@NLO v2.6.1 [44] at LO accuracy
with the NNPDF2.3lo PDF set. The event generator is interfaced with Herwig v7.1.3 [50, 51] to model
the parton shower, hadronization and underlying event. The mass of *X* is varied between 251 GeV and
1000 GeV in the simulation, while its width is set to 10 MeV. The interference with the non-resonant Higgs
boson pair production is neglected.

Single Higgs boson production processes via ggF, VBF, WH, ZH ($qq \to ZH$ and $gg \to ZH$), $t\bar{t}H$, $tH$
($tHqj$ and $tHW$), and $bbH$ are modeled using the same set of MC samples as described in Ref. [52]. For
single Higgs boson production as well as both non-resonant and resonant di-Higgs production, a Higgs boson mass of 125.09 GeV is assumed [53]. The analysis assumes a branching ratio for the Higgs boson decay into two photons of 0.227% and a branching ratio for the Higgs boson decay into two b-quarks of 58.2% [54, 55]. The inclusive cross sections of these processes are normalized to the most precise available theoretical values [54].

The $\gamma\gamma$ plus jets process is simulated with the SHERPA v2.2.4 [56] generator. QCD NLO-accurate matrix elements for up to one parton, and LO-accurate matrix elements for up to three partons are calculated with the Comix [57] and OPENLOOPS 1 [58–60] libraries. These are calculated in the 5-flavor-scheme including b-quarks in the massless approximation and merged with the SHERPA parton shower [61] using the MEPS NLO prescription [62, 63] with a dynamic merging cut [64] of 10 GeV. Within the parton shower, b-quarks are then treated as massive. Finally, events from $t\bar{t}\gamma\gamma$ processes are produced with MADGRAPH5_aMC@NLO in the four-flavor scheme [44]. The simulation samples used in the analysis are listed in Table 1.

<table><caption>Table 1: Summary of single Higgs boson background samples, split by production modes, and continuum background samples. The generator used in the simulation, the PDF set, and tuned parameters (tune) are also provided.</caption><thead><tr><th>Process</th><th>Generator</th><th>PDF set</th><th>Showering</th><th>Tune</th></tr></thead><tbody><tr><td>ggF</td><td>NNLOPS [65–67] [68, 69]</td><td>PDFLHC [42]</td><td>Pythia 8.2 [70]</td><td>AZNLO [71]</td></tr><tr><td>VBF</td><td>POWHEG Box v2 [39, 66, 72–78]</td><td>PDFLHC</td><td>Pythia 8.2</td><td>AZNLO</td></tr><tr><td>WH</td><td>POWHEG Box v2</td><td>PDFLHC</td><td>Pythia 8.2</td><td>AZNLO</td></tr><tr><td>qq → ZH</td><td>POWHEG Box v2</td><td>PDFLHC</td><td>Pythia 8.2</td><td>AZNLO</td></tr><tr><td>gg → ZH</td><td>POWHEG Box v2</td><td>PDFLHC</td><td>Pythia 8.2</td><td>AZNLO</td></tr><tr><td>t<span style="text-decoration: overline;">t</span>H</td><td>POWHEG Box v2 [73–75, 78, 79]</td><td>NNPDF3.0nlo [80]</td><td>Pythia 8.2</td><td>A14 [81]</td></tr><tr><td>bbH</td><td>POWHEG Box v2</td><td>NNPDF3.0nlo</td><td>Pythia 8.2</td><td>A14</td></tr><tr><td>tHqj</td><td>MADGRAPH5_aMC@NLO</td><td>NNPDF3.0nlo</td><td>Pythia 8.2</td><td>A14</td></tr><tr><td>tHW</td><td>MADGRAPH5_aMC@NLO</td><td>NNPDF3.0nlo</td><td>Pythia 8.2</td><td>A14</td></tr><tr><td>γγ+jets</td><td>SHERPA v2.2.4 [56]</td><td>NNPDF3.0nnlo</td><td>SHERPA v2.2.4</td><td>–</td></tr><tr><td>t<span style="text-decoration: overline;">t</span>γγ</td><td>MADGRAPH5_aMC@NLO</td><td>NNPDF2.3lo</td><td>Pythia 8.2</td><td>–</td></tr></tbody></table>

Different pile-up conditions from the same and neighboring bunch crossings are simulated by overlaying the hard-scattering event with inelastic pp events generated by Pythia 8.186 using the NNPDF2.3lo PDF set and the A3 tune [82]. Differences between the simulated and observed distributions of the number of interactions per bunch crossing are corrected by applying pile-up scale factors to simulated events. A full simulation of the ATLAS detector [83] based on GEANT4 [84] is used to reproduce the detector response to single Higgs boson processes. The continuum background and the signal samples are simulated using ATLFASTII [85], a fast simulation of the ATLAS detector response which was shown to be able to accurately simulate di-photon events.
# 4 Object and event selections

## 4.1 Object selection

Photons are reconstructed from topologically-connected clusters [86] of energy deposits in the electromagnetic calorimeter in the region $|η| < 2.37$. The transition region between the barrel and endcap electromagnetic calorimeters, $1.37 < |η| < 1.52$, is excluded. Photon candidates matched to a conversion vertex or a track, which are consistent with originating from photon conversions, are classified as converted photons. Those without a matched conversion vertex or track are classified as unconverted photons.

The calibration of the photon energy is based on a multivariate regression algorithm trained with MC samples, where the input variables are corrected with data-driven techniques. The calibrated energy is finally corrected by applying scale factors derived from $Z \to e^+e^-$ events [87]. The photon direction is reconstructed using the longitudinal segmentation of the calorimeters and constrained to be compatible with the luminous region of the proton beam collisions. Additionally, the conversion vertex position is considered in the case of converted photons.

Events are required to have at least one reconstructed primary vertex, defined as a vertex associated with at least two tracks with transverse momentum ($p_T$) larger than 0.5 GeV.

The primary vertex is selected from the reconstructed collision vertices using a neural-network algorithm [88] based on the extrapolated photon trajectories and the tracks associated to each candidate vertex.

Photon identification is based on the lateral shower profile of the energy deposits in the first and second electromagnetic calorimeter layers and on the energy leakage fraction in the hadronic calorimeter. It reduces the mis-identification of hadronic jets containing large neutral components, primarily $\pi^0$ particles, which decay into a pair of highly collimated photons. 'Tight' identification criteria are applied [87], which are tuned for converted and unconverted photons separately.

To further improve the rejection of mis-identified photons, two isolation variables are defined to quantify the activity around a photon. Calorimeter-based isolation $E_T^{\text{iso}}$ is defined as the sum of the transverse energy of topological clusters within a cone of size $\Delta R = 0.2$ around the photon, correcting for the energy of the photon candidate itself as well as an average expected pile-up contribution. Track-based isolation $p_T^{\text{iso}}$ is defined as the scalar sum of the transverse momenta of all tracks with $p_T > 1$ GeV that originate from the primary vertex, and are within a cone of $\Delta R = 0.2$ around the photon. Isolated photons must have: $E_T^{\text{iso}} < 0.065 \cdot E_T$ and $p_T^{\text{iso}} < 0.05 \cdot E_T$ ('Loose' working point [87]), where $E_T$ is the transverse energy of the photon. For isolated photons and transverse momenta between 30 GeV and 250 GeV, the identification efficiency for unconverted and converted photons ranges from 85% to 99% [87].

Electrons are reconstructed from energy deposits measured in the electromagnetic calorimeter which are matched to ID tracks [87]. They are required to satisfy $|η| < 2.47$, excluding the calorimeter transition region $1.37 < |η| < 1.52$, and have a transverse momentum $p_T > 10$ GeV. Electrons are required to satisfy a 'Medium' identification criterion based on the use of shower shape, track-cluster matching and TRT parameters in a likelihood-based algorithm [87]. Muons are reconstructed from high-quality tracks found in the MS [89]. A matching of these tracks to ID tracks is required in the region $|η| < 2.5$. Muons are required to have $|η| < 2.7$ and $p_T > 10$ GeV, and to satisfy a 'Medium' identification criterion [90]. Both the electrons and muons are matched to the primary vertex via requirements on the tracks' longitudinal and transverse impact parameters, $|z_0|$ and $|d_0|$.
Jets are reconstructed based on particle flow objects built from noise-suppressed positive-energy topological clusters in the calorimeter and reconstructed tracks [91]. The anti-$k_t$ algorithm [92, 93] with a radius parameter of $R = 0.4$ is used. They are required to have $|y| < 4.4$ and $p_T > 25$ GeV. To further suppress jets produced in concurrent pp interactions, each jet within the tracking acceptance of $|η| < 2.4$, and with $p_T < 60$ GeV, is required to satisfy the 'Tight' jet vertex tagger [94] criteria used to identify the jet as originating from the selected primary vertex of the event.

The flavor of jets is determined using a deep-learning neural network, DL1r. The DL1r b-tagging is based on distinctive features of b-hadron decays in terms of the impact parameters of tracks and the displaced vertices reconstructed in the inner detector [95]. The inputs of the DL1r network also include discriminating variables constructed by a recurrent neural network (RNNIP) [96], which exploits the spatial and kinematic correlations between tracks originating from the same b-hadron. This approach is found to improve the performance for jets of high $p_T$ [97] with respect to a previously used multivariate technique [98]. Operating points are defined by a single selection value on the discriminant output distribution and are chosen to provide a specific b-jet efficiency for an inclusive $t\bar{t}$ MC sample. The b-tagging requirements result in an efficiency of 77% for jets containing b-hadrons, and the mis-identification rate is 1/130 (1/4.9) for light-flavor (charm) jets [95]. Scale factors are applied to the simulated events to correct for differences in the b-tagging efficiency between data and simulation. These scale factors are measured as a function of jet $p_T$ using a likelihood-based method in a sample highly enriched with $t\bar{t}$ events [95]. Only jets with $|η| < 2.5$ are considered for flavor-tagging.

The energy of b-tagged jets is corrected for the eventual contribution of muons from semileptonic b-hadron decays. In addition, the undetected energy of the neutrinos and out-of-cone effects are corrected with scale factors derived as a function of the b-jet $p_T$ from a $t\bar{t}$ MC sample. The two corrections together improved the resolution of the invariant mass built with the two jets with the highest b-tagging score ($m_{b\bar{b}}$) by about 20%. The procedure closely follows the one used in Ref. [99], where more details can be found.

An overlap removal procedure is applied to avoid multiple usage of the same detector signals in the same event. Photons are prioritized in this analysis, requiring the removal of electrons, muons and jets within $\Delta R = 0.4$ of a selected photon. Next, jets within $\Delta R = 0.2$ of electrons are removed. In the last step, electrons and muons within $\Delta R = 0.4$ of any jet are removed.

The missing transverse momentum $E_T^{\text{miss}}$ is calculated as the magnitude of the negative vectorial sum of the transverse momenta of all selected and calibrated physics objects of an event that can be matched to the primary vertex. A component called 'soft term' is calculated from the residual tracks that originate from the primary vertex but are not associated with any other object and is added to the $E_T^{\text{miss}}$ calculation [100].

## 4.2 Event selection

Both the non-resonant and resonant searches employ multivariate analysis techniques to select events. Events are selected if they satisfy a common set of preselection requirements; then they are required to fulfill different requirements for the non-resonant and the resonant search.

### 4.2.1 Common preselection

For both the non-resonant and resonant analyses, events are selected using di-photon triggers requiring two reconstructed photon candidates with minimum transverse energies of 35 GeV for the leading and
25 GeV for the sub-leading photon, where leading (sub-leading) refers to the photon candidate with the highest (second-highest) transverse energy [101]. Different photon identification criteria are used: 2015 and 2016 triggers require both photons to fulfill a ‘Loose’ photon identification criteria, while the ‘Medium’ criteria [87] is adopted in 2017–2018 to cope with the increased interaction rate. The trigger efficiency for the non-resonant signal is 82.9% and 69.5% for the resonant signal (using as reference $m_X = 300$ GeV).

On top of the trigger requirements, events are selected if:

* there are at least two photons passing the object selection criteria detailed in Section 4.1;

* the di-photon invariant mass, built with the two leading photons, satisfies $105 \text{ GeV} < m_{\gamma\gamma} < 160 \text{ GeV}$;

* the leading (sub-leading) photon $p_T$ is larger than 35% (25%) of the mass of the di-photon system;

* there are exactly two b-tagged jets;

* no electrons or muons are present;

* fewer than six central ($|\eta| < 2.5$) jets are required, to help rejecting $t\bar{t}H$ events where the top quarks decay hadronically.

In order to remain orthogonal to the ATLAS search for $HH \to b\bar{b}b\bar{b}$ [102], any event with more than two b-jets using the 77% efficient working point is rejected. Multivariate techniques are then used to target the non-resonant ggF production mode or the resonant production mode of Higgs boson pairs. For both the non-resonant and resonant analyses, the resulting boosted decision tree (BDT)-based categories are required to have at least 9 expected continuum background events in the range $105 \text{ GeV} < m_{\gamma\gamma} < 160 \text{ GeV}$ excluding $120 \text{ GeV} < m_{\gamma\gamma} < 130 \text{ GeV}$, which guarantees that the final selection on the data sideband has sufficient events to perform an $m_{\gamma\gamma}$ fit.

The $m^*_{b\bar{b}\gamma\gamma}$ variable, defined as $m^*_{b\bar{b}\gamma\gamma} = m_{b\bar{b}\gamma\gamma} - m_{b\bar{b}} - m_{\gamma\gamma} + 250 \text{ GeV}$ (where 250 GeV is about twice the Higgs boson mass value), is used to implement selection criteria for both the non-resonant and resonant analyses. Compared with $m_{b\bar{b}\gamma\gamma}$, the $m^*_{b\bar{b}\gamma\gamma}$ variable improves the four-object resolution, in particular for resonant signals decaying into a pair of Higgs bosons, as shown in Figure 4 thanks to cancellation of detector resolution effects.

### 4.2.2 Non-resonant selection

Following the preselection, events are divided into two regions using the value of the $m^*_{b\bar{b}\gamma\gamma}$ variable. A high mass region, with $m^*_{b\bar{b}\gamma\gamma} > 350 \text{ GeV}$, targets the SM signal ($\kappa_\lambda = 1$), while a low mass region, with $m^*_{b\bar{b}\gamma\gamma} < 350 \text{ GeV}$, is used to retain sensitivity for BSM signals ($\kappa_\lambda = 10$). The dependence of $m^*_{b\bar{b}\gamma\gamma}$ on $\kappa_\lambda$ can be seen in Figure 5.

In each mass region, a dedicated BDT is trained using XGBoost [103] to discriminate between a benchmark $HH$ signal and a combination of $\gamma\gamma$, $t\bar{t}H$, $ggH$, and $ZH$ simulated backgrounds. In the high mass region, the SM $HH$ sample is used as signal, while in the low mass region, the $\kappa_\lambda = 10$ sample is used as signal.

The BDT input variables are summarized in Table 2. Identical variable sets are used for high mass and low mass categories.

The BDT combines several input variables that exploit the different kinematic properties of signal and background events, as well as the b-tagging information. Observables based on the kinematic properties of
<figure><img src="image_5.png" /><figcaption>Figure 4: Reconstructed four-body mass for $m_X = 300$ GeV and $m_X = 500$ GeV resonant signal benchmarks and for the $\gamma\gamma$+jets background. Dashed lines represent the distribution of $m_{b\bar{b}\gamma\gamma}$ while solid lines represent the distribution of $m^*_{b\bar{b}\gamma\gamma}$, defined in Section 4.2.1. Distributions are normalized to unit area.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: The $m^*_{b\bar{b}\gamma\gamma}$ distributions after the common preselection for (a) non-resonant ggF $HH$ and (b) VBF $HH$ signals with several $\kappa_\lambda$ values. $m^*_{b\bar{b}\gamma\gamma} = 350$ GeV is chosen as the separating boundary between categories targeting the SM and BSM $\kappa_\lambda$ signals.</figcaption></figure>

the reconstructed photons, such as the leading and sub-leading photon's angular information, the transverse momentum over the invariant mass of the di-photon system are combined with jet-based information. The 'single topness' variable ($\chi_{W_t}$) is also used. It is defined in Equation (1).

$$ \chi_{W_t} = \min \sqrt{\left(\frac{m_{j_1 j_2} - m_W}{m_W}\right)^2 + \left(\frac{m_{j_1 j_2 j_3} - m_t}{m_t}\right)^2}, \qquad (1) $$

where the minimum is taken over all combinations of three jets in the event (with no requirements on b-tagging status), $m_W = 80$ GeV, and $m_t = 173$ GeV. Among the input variables, $m_{b\bar{b}}$ and $H_T$ show the highest discriminating power against the $\gamma\gamma$ plus jets continuum background. Particular care has been taken to avoid the BDT event selection leading to biases in the $m_{\gamma\gamma}$ background distribution. The BDT score
distributions of the low mass and high mass regions are shown in Figure 6 for events passing the common preselection. In each mass region, two categories are defined based on the BDT score. The boundaries of the categories are chosen by maximizing the combined number counting significance [104] using signal and background yields in the mass window $120 \text{ GeV} < m_{\gamma\gamma} < 130 \text{ GeV}$ in the chosen categories. The four resulting BDT categories are defined in Table 3.

<table><caption>Table 2: Variables used in the BDT for the non-resonant analysis. The b-tag status identifies the highest fixed b-tag working point (60%, 70%, 77%) that the jet passes. All vectors in the event are rotated so that the leading photon $\phi$ is equal to zero.</caption><thead><tr><th>Variable</th><th>Definition</th></tr></thead><tbody><tr><td colspan="2">Photon-related kinematic variables</td></tr><tr><td>p<sub>T</sub>/m<sub>γγ</sub></td><td>Transverse momentum of the two photons scaled by their invariant mass m<sub>γγ</sub></td></tr><tr><td>η and φ</td><td>Pseudo-rapidity and azimuthal angle of the leading and sub-leading photon</td></tr><tr><td colspan="2">Jet-related kinematic variables</td></tr><tr><td>b-tag status</td><td>Highest fixed b-tag working point that the jet passes</td></tr><tr><td>p<sub>T</sub>, η and φ</td><td>Transverse momentum, pseudo-rapidity and azimuthal angle of the two jets with the highest b-tagging score</td></tr><tr><td>p<sub>T</sub><sup>b<span style="text-decoration: overline">b</span></sup>, η<sub>b<span style="text-decoration: overline">b</span></sub> and φ<sub>b<span style="text-decoration: overline">b</span></sub></td><td>Transverse momentum, pseudo-rapidity and azimuthal angle of b-tagged jets system</td></tr><tr><td>m<sub>b<span style="text-decoration: overline">b</span></sub></td><td>Invariant mass built with the two jets with the highest b-tagging score</td></tr><tr><td>H<sub>T</sub></td><td>Scalar sum of the p<sub>T</sub> of the jets in the event</td></tr><tr><td>Single topness</td><td>For the definition, see Eq. (1)</td></tr><tr><td colspan="2">Missing transverse momentum-related variables</td></tr><tr><td>E<sub>T</sub><sup>miss</sup> and φ<sup>miss</sup></td><td>Missing transverse momentum and its azimuthal angle</td></tr></tbody></table>

<table><caption>Table 3: Definition of the categories used in the HH non-resonant search. Before entering the BDT-based categories, events are required to satisfy the common preselection.</caption><thead><tr><th>Category</th><th>Selection criteria</th></tr></thead><tbody><tr><td>High mass BDT tight</td><td>m<sup>*</sup><sub>b<span style="text-decoration: overline">b</span>γγ</sub> ≥ 350 GeV, BDT score ∈ [0.967, 1]</td></tr><tr><td>High mass BDT loose</td><td>m<sup>*</sup><sub>b<span style="text-decoration: overline">b</span>γγ</sub> ≥ 350 GeV, BDT score ∈ [0.857, 0.967]</td></tr><tr><td>Low mass BDT tight</td><td>m<sup>*</sup><sub>b<span style="text-decoration: overline">b</span>γγ</sub> < 350 GeV, BDT score ∈ [0.966, 1]</td></tr><tr><td>Low mass BDT loose</td><td>m<sup>*</sup><sub>b<span style="text-decoration: overline">b</span>γγ</sub> < 350 GeV, BDT score ∈ [0.881, 0.966]</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 6: The BDT distribution of the di-Higgs ggF signal for two different values of $k_\lambda$ and the main backgrounds in the (a) low and (b) high mass region. Distributions are normalized to unit area. The dotted lines denote the category boundaries. Events with a BDT score below 0.881 in the low mass region or below 0.857 in the high mass region are discarded.</figcaption></figure>

### 4.2.3 Resonant selection

The resonant analysis uses a multivariate analysis based on a BDT technique. A potential limitation of a BDT-based selection is the low background statistics for higher resonance masses. To overcome this limitation, a single BDT is trained for all resonance masses with the signal reweighted event-by-event to match the $m^*_{b\bar{b}\gamma\gamma}$ distribution of the background events, such that the training is independent of the resonant signal mass hypothesis.

Using the TMVA toolkit [105], two BDTs are trained to better separate the signal from backgrounds of different nature: the $\gamma\gamma$ plus the $t\bar{t}\gamma\gamma$ backgrounds (BDT$_{\gamma\gamma}$) and the single Higgs boson background (BDT$_{SingleH}$), where ZH and $t\bar{t}H$ production modes are the dominant resonant backgrounds. A complete list of the variables used for the BDT training is detailed in Table 4. The $E_T^{\text{miss}}$ information is also used in the training as it is useful to reject the single Higgs boson ($t\bar{t}H$ in particular) and the $t\bar{t}\gamma\gamma$ backgrounds.

The combined BDT score of an event is obtained by combining the two BDT scores in quadrature, as shown in Equation (2):

$$ \text{BDT}_{\text{tot}} = \frac{1}{\sqrt{C_1^2 + C_2^2}} \sqrt{C_1^2 \left( \frac{\text{BDT}_{\gamma\gamma} + 1}{2} \right)^2 + C_2^2 \left( \frac{\text{BDT}_{\text{SingleH}} + 1}{2} \right)^2} \quad (2) $$

The coefficients $C_1, C_2$ ($C_2 = 1 - C_1$) and the BDT$_{tot}$ take values in the range [0,1]. Only events passing a minimum requirement on the value of the BDT$_{tot}$ are considered in the analysis. The values of the $C_1$ and $C_2$ as well as the BDT cut value are scanned twice in order to maximize the significance. This two-stage optimization procedure finds first the maximum significance that can be achieved for each resonance mass point independently, leading to different coefficients and BDT$_{tot}$ requirement value for each mass parameter. A second scan is done to select all coefficients providing a significance within 5% from the maximum value, for each of the resonance mass value. From those possible combinations a common $C_1$ coefficient
<table><caption>Table 4: Variables used in the BDT for the resonant analysis. For variables depending on b-tagged jets, only jets b-tagged using the 77% working point are considered as described in Section 4.1.</caption><thead><tr><th>Variable</th><th>Definition</th></tr></thead><tbody><tr><td colspan="2"><strong>Photon-related kinematic variables</strong></td></tr><tr><td>p<sup>γγ</sup><sub>T</sub>, y<sup>γγ</sup></td><td>Transverse momentum and rapidity of the di-photon system</td></tr><tr><td>Δφ<sub>γγ</sub> and ΔR<sub>γγ</sub></td><td>Azimuthal angular distance and ΔR between the two photons</td></tr><tr><td colspan="2"><strong>Jet-related kinematic variables</strong></td></tr><tr><td>m<sub>b<span style="text-decoration: overline">b</span></sub>, p<sup>b<span style="text-decoration: overline">b</span></sup><sub>T</sub> and y<sub>b<span style="text-decoration: overline">b</span></sub></td><td>Invariant mass, transverse momentum and rapidity of the b-tagged jets system</td></tr><tr><td>Δφ<sub>b<span style="text-decoration: overline">b</span></sub> and ΔR<sub>b<span style="text-decoration: overline">b</span></sub></td><td>Azimuthal angular distance and ΔR between the two b-tagged jets</td></tr><tr><td>N<sub>jets</sub> and N<sub>b-jets</sub></td><td>Number of jets and number of b-tagged jets</td></tr><tr><td>H<sub>T</sub></td><td>Scalar sum of the p<sub>T</sub> of the jets in the event</td></tr><tr><td colspan="2"><strong>Photons and jets-related kinematic variables</strong></td></tr><tr><td>m<sub>b<span style="text-decoration: overline">b</span>γγ</sub></td><td>Invariant mass built with the di-photon and b-tagged jets system</td></tr><tr><td>Δy<sub>γγ,b<span style="text-decoration: overline">b</span></sub>, Δφ<sub>γγ,b<span style="text-decoration: overline">b</span></sub> and ΔR<sub>γγ,b<span style="text-decoration: overline">b</span></sub></td><td>Distance in rapidity, azimuthal angle and ΔR between the di-photon and the b-tagged jets system</td></tr></tbody></table>

$(C_1 = 0.65)$ is searched across all the resonances so that the selection will have common coefficients for all resonance mass points, but different BDT<sub>tot</sub> values. For each of the mass hypotheses of the resonance, a requirement is set on the $m^*_{b\bar{b}\gamma\gamma}$ value to select events within $\pm 2\sigma$ of the expected mean value for signal events (where $\sigma$ is the standard deviation parameter of a fit of the $m^*_{b\bar{b}\gamma\gamma}$ distribution with a Crystal Ball function). In the case of the 900 GeV and 1000 GeV mass hypotheses, the requirement is relaxed to $\pm 4\sigma$ to increase the number of events used for the signal extraction. The BDT<sub>tot</sub> distributions are shown in Figure 7.

## 4.3 Data/Predictions comparison after selection

The analysis selection described in Section 4.2 requires two tight photons and this region is mainly composed of $\gamma\gamma$, $\gamma$-jet, and di-jet events, where either one or two jets are mis-identified as photons. The fraction of each component can be estimated using a data-driven approach [106] using the photon identification and isolation distributions from genuine and mis-identified photons. After the common preselection, 85±3% of sideband events consist of genuine di-photon events, with the remaining 15±4% consisting of $\gamma$-jet events and a negligible amount of di-jet events. The uncertainties on the above fractions are calculated considering both statistics and systematic uncertainties, where the systematic uncertainty is estimated using different photon identification criteria. In the BDT-based categories, the fraction of
<figure><img src="image_8.png" /><figcaption>Figure 7: The BDT score for the benchmark signals ((a) $m_X$ = 300 GeV and (b) $m_X$ = 500 GeV) and for the main backgrounds. Distributions are normalized to unit area. The dotted lines denote the event selection thresholds. Events with a BDT score below 0.85 for $m_X$ = 300 GeV or below 0.75 for $m_X$ = 500 GeV are discarded.</figcaption></figure>

genuine di-photon events increases but the method suffers from low statistics for both the non-resonant and resonant cases.

Figure 8 shows the agreement between data and background prediction for the $m_{γγ}$ and $m_{b\bar{b}γγ}^*$ distributions, after the common preselection. The continuum background is scaled by the γγ, γ-jet, and di-jet fractions and normalized to the data sideband. The γγ+jets continuum background is further divided according to the flavors of the two jets (for example *bb* or other jets). This decomposition is taken directly from the proportions predicted by the SHERPA event generator, as described in Section 3, and it is shown for illustration. Figures 9 and 10 show the $m_{γγ}$ distribution after the non-resonant and resonant BDT categorization and for two benchmark mass points $m_X$ = 300 GeV and $m_X$ = 500 GeV. The figures are illustrative to show the signal and background composition. These studies are not used to determine the background in the analysis workflow.

# 5 Signal and background parametrization

The signal and backgrounds are extracted by fitting analytic functions to the di-photon invariant mass distribution in the range $105 \text{ GeV} < m_{γγ} < 160 \text{ GeV}$ in both the resonant and non-resonant *HH* searches.

## 5.1 Signal parameterization

For the di-Higgs signal and the single Higgs boson production background processes, the parameterized forms are determined through fits to simulated samples and the expected normalizations are obtained from their theoretical cross sections multiplied by the product of the acceptance times efficiency from the simulation. The di-photon invariant mass distribution shapes are modeled with a double-sided Crystal
<figure><img src="image_9.png" /><figcaption>Figure 8: Distributions of (a) $m_{\gamma\gamma}$ and (b) $m_{b\bar{b}\gamma\gamma}^*$ for events passing the common preselection criteria. The continuum background is scaled by the $\gamma\gamma$, $\gamma$-jet or jet-$\gamma$, and di-jet fractions and normalized to the data sideband.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Distributions of $m_{\gamma\gamma}$ in all signal categories for the non-resonant $\text{HH}$ search: (a) high mass BDT tight, (b) high mass BDT loose, (c) low mass BDT tight, (d) low mass BDT loose. The continuum background is scaled by the $\gamma\gamma$, $\gamma$-jet, and di-jet fractions and normalized to the data sideband.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Distributions of $m_{γγ}$ for the selections used for the resonance mass points (a) $m_X = 300$ GeV and (b) $m_X = 500$ GeV for the resonant search. The non-resonant background is scaled by the $γγ$, $γ$-jet, and di-jet fractions and normalized to the data sideband. The scalar resonance signal is scaled to an arbitrary cross section value.</figcaption></figure>

Ball function [88, 107], which is characterized by a Gaussian core and asymmetric power law tails. This function allows the modeling of situations in which non-Gaussian tails can arise from experimental effects, such as photon energy mis-measurements.

The shape parameters are determined by fitting the di-photon mass distribution in simulation for each category. The width of the fitted function is largely insensitive to the specific signal processes considered in the analysis, with maximum variations of approximately 10%. For the non-resonant search, the parameterized form of $m_{γγ}$ is obtained from the simulation of the ggF and VBF $HH$ processes with $\kappa_{\lambda} = 1$, described in Section 3. No significant dependence of the functional form with $\kappa_{\lambda}$ was found. For the resonant search, the functional form is obtained from the simulation of the heavy resonance signals. Table 5 shows the effective resolution (the smallest mass window that contains 68% of the signal events, termed $\sigma_{68}$) of the functional form fit to the $m_{γγ}$ distribution for simulated Higgs boson pair events for the non-resonant categories and for two different mass hypotheses for a heavy resonant signal. For both searches, the chosen functional forms are found to model both the single Higgs and di-Higgs boson events well. As no statistically significant bias is observed in injection tests between the input and fitted signals, the same parameterized functions are used.

<table><caption>Table 5: Effective resolution of the $m_{γγ}$ invariant mass spectrum (the smallest mass window that contains 68% of signal, termed $\sigma_{68}$) and corresponding statistical uncertainty are given for the non-resonant search categories and for the two benchmark scalar resonant signals.</caption><thead><tr><th>Category</th><th>σ<sub>68</sub> [GeV]</th></tr></thead><tbody><tr><td>High mass BDT tight</td><td>1.46 ± 0.01</td></tr><tr><td>High mass BDT loose</td><td>1.61 ± 0.02</td></tr><tr><td>Low mass BDT tight</td><td>1.72 ± 0.06</td></tr><tr><td>Low mass BDT loose</td><td>1.81 ± 0.03</td></tr><tr><td>Resonant m<sub>X</sub> = 300 GeV</td><td>1.96 ± 0.02</td></tr><tr><td>Resonant m<sub>X</sub> = 500 GeV</td><td>1.60 ± 0.01</td></tr></tbody></table>
## 5.2 Background parameterization

The continuum di-photon background is modeled using a functional form chosen by fitting a high statistics MC background template. Given the high $\gamma\gamma$ purity quoted in Section 4.3, the background template is constructed from the $\gamma\gamma$+jets simulation which is normalized to the data sideband in the mass windows of 105–120 GeV and 130–160 GeV in the $m_{\gamma\gamma}$ spectrum. The difference in shape between simulated events and those measured in the sidebands is found to have a negligible impact on the background parametrization.

The potential bias associated with the choice of a specific analytic function to model the continuum background is estimated for each category as the signal event yield extracted from a signal-plus-background fit to the background-only di-photon invariant mass distribution in the range $105 \text{ GeV} < m_{\gamma\gamma} < 160 \text{ GeV}$, as prescribed in Refs. [88, 108]. This bias is also termed ‘spurious signal’. The number of fitted signal events is computed for Higgs boson masses varying in intervals of 1 GeV from 121 GeV to 129 GeV. The bias is taken as the largest number of fitted signal events in this 8 GeV mass window. Of the different analytic functions that are tested, the one with the smallest number of parameters is chosen among the functions with a spurious signal smaller than 20% of the data statistical uncertainty plus two times of the MC statistical uncertainty. In each category of the non-resonant and in all the analysis regions defined in the resonant HH search, an exponential function $\exp(a \cdot m_{\gamma\gamma})$ is found to be the best choice since it has the minimal number of degrees of freedom and yields a consistently small bias. Wald tests [109] on data show that the data do not prefer a higher degree functional form with respect to the exponential form.

# 6 Systematic uncertainties

The sensitivity of the analysis is limited by the statistical precision. The assessment of the systematic uncertainties is described below. The uncertainty in the integrated luminosity of the full Run 2 data set is 1.7% [37], obtained using the LUCID-2 detector [36] for the primary luminosity measurements.

The main background processes of the analysis are estimated from data and are subject to uncertainties related to the potential bias arising from the chosen background model, as described in section 5.2. The background functional form bias is assessed as an additional uncertainty in the total number of signal events in each category. For the single Higgs boson and di-Higgs boson production processes, both of which are estimated using simulation, experimental and theoretical systematic uncertainties are considered through the full analysis procedure, as described in the following.

The efficiency of the di-photon trigger used to select events is evaluated in simulation using a trigger matching technique and in data using a bootstrap method [101]. In the di-photon invariant mass range of $105 \text{ GeV} < m_{\gamma\gamma} < 160 \text{ GeV}$, the trigger efficiency uncertainty affects the acceptance by 1% in each category. The uncertainty in the vertex selection efficiency is assessed by comparing the efficiency of finding photon-pointing vertices in $Z \to e^+e^-$ events in data with that in simulation [110]. The resulting uncertainty is found to have a negligible effect on the signal selection efficiency.

The systematic uncertainties due to the photon identification and isolation efficiencies are estimated following the prescriptions in Ref. [87]. They are evaluated by varying the correction factors of photon selection efficiencies in simulation by the corresponding uncertainties and affect the di-photon selection efficiency. The experimental uncertainties in photon scale and resolution are obtained from Ref. [87].
The jet energy scale and resolution uncertainties affect the $m_{b\bar{b}}$ distribution, while flavor-tagging uncertainties affect the acceptance of the analysis categories. The experimental uncertainties in jet energy scale and resolution are propagated to the $E_T^{\text{miss}}$ calculation. In addition, the uncertainties in the scale and resolution of the $E_T^{\text{miss}}$ soft term are evaluated by using the method described in Ref. [100]. The flavor-tagging uncertainties for b- and c-jets are estimated in $t\bar{t}$ events [95, 111], while the mis-identification uncertainty of light-flavor jets is determined using di-jet events [112]. Additional uncertainties related to the b-jet momentum correction related to the eventual presence of muons and neutrinos are found to be negligible.

For the single Higgs boson and SM HH production, the effects of theoretical scale uncertainties due to missing higher-order corrections on the production rates are estimated by varying the factorization and renormalization scales up and down from their nominal values by a factor of two, recalculating the cross section in each case, and taking the largest deviation from the nominal cross section as the uncertainty. The uncertainties in the cross sections, including the PDF+$\alpha_s$ uncertainties, and the uncertainties on the $H \to \gamma\gamma$ and $H \to b\bar{b}$ branching fractions are taken from Ref. [13, 54]. The uncertainty on the value of the Higgs boson mass is considered [53]. A 100% uncertainty is assigned to the single Higgs boson ggF and VBF production modes and to the production in association with a W boson. This is motivated by studies of heavy-flavor production in association with top-quark pairs [113, 114] and of W bosons in association with b-jets [115]. No additional heavy-flavor uncertainty is assigned to single Higgs boson $t\bar{t}H$ and ZH production modes, where the dominant heavy-flavor production is already accounted for in the LO process. In addition, for the non-resonant HH production processes, a systematic uncertainty is assigned to the $\kappa_\lambda$ reweighting, as discussed in Section 3.

In the resonant search, uncertainties arising from scale uncertainties are neglected. In this search the SM HH production processes are considered as background. For all samples, uncertainties related to the PDF choice and differences between alternative models of parton showering and hadronization are considered.

# 7 Results

The statistical framework used to derive the results for both the non-resonant and resonant searches are described in the following.

## 7.1 Statistical framework

For both the non-resonant and resonant searches, the results of the analysis are obtained from a maximum likelihood fit of the $m_{\gamma\gamma}$ distribution in the range $105 \text{ GeV} < m_{\gamma\gamma} < 160 \text{ GeV}$, performed simultaneously over all relevant categories described in Section 4.2. The likelihood function is defined in Eq. (3):

$$ \mathcal{L} = \prod_c \left( \text{Pois}(n_c | N_c(\boldsymbol{\theta})) \cdot \prod_{i=1}^{n_c} f_c(m_{\gamma\gamma}^i, \boldsymbol{\theta}) \cdot G(\boldsymbol{\theta}) \right), \quad (3) $$

where for each event $i$ in a category $c$, $n_c$ is the observed number of events, $N_c$ is the expected number of events, $f_c$ is the value of the probability density function (pdf), $\boldsymbol{\theta}$ are nuisance parameters, and $G(\boldsymbol{\theta})$ are constraint pdfs on the nuisance parameters.
The expected number of events $N_c$, defined in Eq. (4), is the sum of the expected yields from di-Higgs boson production processes, single Higgs boson, the non-resonant background, as well as the spurious signal uncertainty.

$$N_c(\boldsymbol{\theta}) = \mu \cdot N_{HH,c}(\boldsymbol{\theta}_{HH}^{\text{yield}}) + N_{bkg,c}^{\text{res}}(\boldsymbol{\theta}_{res}^{\text{yield}}) + N_{SS,c} \cdot \boldsymbol{\theta}^{\text{SS,c}} + N_{bkg,c}^{\text{non-res}}, \quad (4)$$

where $\mu$ is the signal strength, $\boldsymbol{\theta}^{\text{SS,c}}$ represent the nuisance parameters associated to the background function bias and $\boldsymbol{\theta}^{\text{yield}}$ represent the nuisance parameters affecting the events yield, as detailed in Section 6. Correlation of the nuisance parameters across different signal and background components, as well as categories, is taken into account. In the case of the non-resonant search $N_{bkg,c}^{\text{res}} = N_{H,c}$, while in the case of the resonant analysis $N_{bkg,c}^{\text{res}} = N_{H,c} + N_{SMHH,c}$.

The probability density function $f_c$ representing the shape information. The sum of the double-sided Crystal Ball functions modeling the HH production processes, single Higgs boson, and the spurious signal, and of the analytical function modeling the non-resonant background as described in Section 5.2 is shown in Eq. (5):

$$f_c(m_{\gamma\gamma}, \boldsymbol{\theta}) = [\mu \cdot N_{HH,c}(\boldsymbol{\theta}_{HH}^{\text{yield}}) \cdot f_{HH,c}(m_{\gamma\gamma}, \boldsymbol{\theta}_{HH}^{\text{shape}}) + N_{bkg,c}^{\text{res}}(\boldsymbol{\theta}_{res}^{\text{yield}}) \cdot f_{bkg,c}^{\text{res}}(m_{\gamma\gamma}, \boldsymbol{\theta}_{res}^{\text{shape}}) + N_{SS,c} \cdot \boldsymbol{\theta}_{HH}^{\text{SS,c}} \cdot f_{HH,c}(m_{\gamma\gamma}, \boldsymbol{\theta}_{HH}^{\text{shape}}) + N_{bkg,c}^{\text{non-res}} \cdot f_{bkg,c}^{\text{non-res}}(m_{\gamma\gamma}, \boldsymbol{\theta}_{non-res}^{\text{shape}})]/N_c(\boldsymbol{\theta}_{non-res}^{\text{yield}}), \quad (5)$$

where $\boldsymbol{\theta}^{\text{shape}}$ represents nuisance parameters related to the shape variations of the functional forms. When a nuisance parameter is related to shape and yield variations at the same time, the two effects are correlated.

The nominal yields of the resonant background processes are set to values from simulation. The signal strength, non-resonant background shape parameters, and the nuisance parameters representing the systematic uncertainties are free parameters in the fit. The measurement of the parameters of interest is carried out using a statistical test based on the profile likelihood ratio [104].

In the absence of signal, exclusion limits are set on Higgs boson pair production in the $b\bar{b}\gamma\gamma$ final state. The limits for both non-resonant and resonant production are calculated using the CL$_S$ method [116], with the profile-likelihood-ratio-based test statistic $\tilde{q}_{\mu}$. The asymptotic approximation [104] is used for the test-statistic distribution.

## 7.2 Non-resonant search results

Figure 11 shows the background fits to the data. No significant excess over the SM background expectations is found, as summarized in Table 6. Therefore, limits at 95% CL are set based on the profile likelihood ratio approach. In the absence of signal, the statistical analysis sets a 95% CL upper limit on the non-resonant HH production cross section at 130 fb, while 180 fb is expected. An observed (expected) upper limit at 95% CL on the signal strength of 4.1 (5.5) times the SM prediction is obtained. For the upper limits on the cross section, all theoretical uncertainties are included, except those related to the signal cross section itself, while constraints on the signal strength are computed including uncertainties in the predicted signal cross section. The dominant systematic uncertainties are listed in Table 8. The observed (expected) constraints on $\kappa_{\lambda}$ are [-1.5,6.7] ([-2.4, 7.7]) at 95% CL, as shown in Figure 12. The single Higgs boson production cross sections and Higgs boson decay branching ratios are assumed to have SM values in this paper. The coupling strength between Higgs boson and other particles are set to their SM values [52]. The constraints
on $\kappa_{\lambda}$ are obtained over an expected hypothesis excluding $pp \to HH$ production. The inclusion of the VBF $HH$ production mode improves the constraints by about 5% with respect to an alternative fit considering only the ggF production mode.

<table><caption>Table 6: Expected and observed numbers of events in the categories of the non-resonant search. An additional requirement of 120 GeV < $m_{\gamma\gamma}$ < 130 GeV is applied. The uncertainties on the continuum background are those arising from the fitting procedure. The uncertainties on the single Higgs boson and Higgs boson pair productions are from MC statistical error.</caption><thead><tr><th></th><th>High mass<br>BDT tight</th><th>High mass<br>BDT loose</th><th>Low mass<br>BDT tight</th><th>Low mass<br>BDT loose</th></tr></thead><tbody><tr><td>Continuum background</td><td>4.9 ± 1.1</td><td>9.5 ± 1.5</td><td>3.7 ± 1.0</td><td>24.9 ± 2.5</td></tr><tr><td>Single Higgs boson background</td><td>0.670 ± 0.032</td><td>1.57 ± 0.04</td><td>0.220 ± 0.016</td><td>1.39 ± 0.04</td></tr><tr><td>ggF</td><td>0.261 ± 0.028</td><td>0.44 ± 0.04</td><td>0.063 ± 0.014</td><td>0.274 ± 0.030</td></tr><tr><td>t<span style="text-decoration: overline;">t</span>H</td><td>0.1929 ± 0.0045</td><td>0.491 ± 0.007</td><td>0.1074 ± 0.0033</td><td>0.742 ± 0.009</td></tr><tr><td>ZH</td><td>0.142 ± 0.005</td><td>0.486 ± 0.010</td><td>0.04019 ± 0.0027</td><td>0.269 ± 0.007</td></tr><tr><td>Rest</td><td>0.074 ± 0.012</td><td>0.155 ± 0.020</td><td>0.008 ± 0.006</td><td>0.109 ± 0.016</td></tr><tr><td>SM HH signal</td><td>0.8753 ± 0.0032</td><td>0.3680 ± 0.0020</td><td>(49.4 ± 0.7) · 10<sup>-3</sup></td><td>(78.7 ± 0.9) · 10<sup>-3</sup></td></tr><tr><td>ggF</td><td>0.8626 ± 0.0032</td><td>0.3518 ± 0.0020</td><td>(46.1 ± 0.7) · 10<sup>-3</sup></td><td>(71.8 ± 0.9) · 10<sup>-3</sup></td></tr><tr><td>VBF</td><td>0.01266 ± 0.00016</td><td>0.01618 ± 0.00018</td><td>(3.22 ± 0.08) · 10<sup>-3</sup></td><td>(6.923 ± 0.011) · 10<sup>-3</sup></td></tr><tr><td>Alternative HH($\kappa_{\lambda}$ = 10) signal</td><td>6.36 ± 0.05</td><td>3.691 ± 0.038</td><td>4.65 ± 0.04</td><td>8.64 ± 0.06</td></tr><tr><td>Data</td><td>2</td><td>17</td><td>5</td><td>14</td></tr></tbody></table>

## 7.3 Resonant search results

Figure 13 shows the fit to the data of the resonant search for two benchmark values of the mass $m_X$ of a hypothetical scalar particle. No significant excess over the SM background expectations is found, as shown in Table 7. Figure 14 shows the observed and expected upper limits at 95% CL on the production cross section of a narrow width scalar resonance. The observed (expected) upper limits vary between 610–47 fb (360–43 fb) in the range $251 \text{ GeV} \le m_X \le 1000 \text{ GeV}$.

The dominant systematic uncertainties are listed in Table 8. The main uncertainties are related to the choice of the functional form of the continuum background and to the photon energy scale and resolution.
<figure><img src="image_12.png" /><figcaption>Figure 11: Data are compared to the background-only fit for the four categories of the non-resonant search. Both the continuum background and the background from single Higgs boson production are considered.</figcaption></figure>

<table><caption>Table 7: Expected and observed numbers of events of the resonant *HH* search. An additional requirement of 120 GeV < *m*<sub>γγ</sub> < 130 GeV is applied. The event numbers quoted for the scalar resonance signal assume an arbitrary total production cross section *σ*(*pp* → *X* → *HH*) equal to the observed exclusion limits of Figure 14. The uncertainties on the continuum background are those arising from the fitting procedure. The uncertainties on the single Higgs boson, Higgs boson pair and scalar resonance production are from the MC statistical error.</caption><thead><tr><th></th><th>m<sub>X</sub> = 300 GeV</th><th>m<sub>X</sub> = 500 GeV</th></tr></thead><tbody><tr><td>Continuum background</td><td>5.6 ± 2.4</td><td>3.5 ± 2.0</td></tr><tr><td>Single Higgs boson background</td><td>0.339 ± 0.009</td><td>0.398 ± 0.010</td></tr><tr><td>SM <i>HH</i> background</td><td>(20.6 ± 0.5) · 10<sup>-3</sup></td><td>0.1932 ± 0.0015</td></tr><tr><td><i>X</i> → <i>HH</i> signal</td><td>5.771 ± 0.031</td><td>5.950 ± 0.026</td></tr><tr><td>Data</td><td>6</td><td>4</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Figure 12: Observed and expected limits at 95% CL on the cross section of non-resonant Higgs boson pair production as a function of the Higgs boson self-coupling modifier $\kappa_{\lambda} = \lambda_{HHH}/\lambda_{HHH}^{SM}$. The constraints on $\kappa_{\lambda}$ are obtained over an expected hypothesis excluding $pp \to HH$ production. The $\pm 1\sigma$ and $\pm 2\sigma$ variations about the expected limit due to statistical and systematic uncertainties are also shown. The theory prediction curve represents the scenario where all parameters and couplings are set to their SM values except for $\kappa_{\lambda}$. The uncertainty band of the theory prediction curve shows the cross section uncertainty.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Data are compared to the background-only fit for the resonant search for the (a) $m_X = 300$ GeV and (b) $m_X = 500$ GeV mass hypotheses. The continuum background, as well as the background from single Higgs boson production and from the SM $HH$ production are considered.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: Observed and expected limits at 95% CL on the production cross section of a narrow width scalar resonance X as a function of the mass $m_X$ of the hypothetical scalar particle. The black solid line represents the observed upper limits. The dashed line represents the expected upper limits. The $\pm 1\sigma$ and $\pm 2\sigma$ variations about the expected limit due to statistical and systematic uncertainties are also shown.</figcaption></figure>

<table><caption>Table 8: Breakdown of the dominant systematic uncertainties. The impact of the uncertainties is defined according to the statistical analysis described in Section 7. It corresponds to the variation on the upper limit on the signal strength when re-evaluating the profile likelihood ratio after fixing the nuisance parameter in question to its best-fit value increased or decreased by one standard deviation, while all remaining nuisance parameters remain free to float. The impact is shown in %. Only systematic uncertainties with an impact of at least 0.5% are shown. Uncertainties of Norm. + Shape type have effects on both the normalization and the parameters of the functional form, the rest of uncertainties affects only the yields.</caption><thead><tr><th rowspan="2">Source</th><th rowspan="2">Type</th><th colspan="2">Relative impact of the systematic uncertainties in %</th></tr><tr><th>Non-resonant analysis<br>HH</th><th>Resonant analysis<br>$m_X$ = 300 GeV</th></tr></thead><tbody><tr><td colspan="4">Experimental</td></tr><tr><td>Photon energy scale</td><td>Norm. + Shape</td><td>5.2</td><td>2.7</td></tr><tr><td>Photon energy resolution</td><td>Norm. + Shape</td><td>1.8</td><td>1.6</td></tr><tr><td>Flavor tagging</td><td>Normalization</td><td>0.5</td><td>&lt; 0.5</td></tr><tr><td colspan="4">Theoretical</td></tr><tr><td>Heavy flavor content</td><td>Normalization</td><td>1.5</td><td>&lt; 0.5</td></tr><tr><td>Higgs boson mass</td><td>Norm. + Shape</td><td>1.8</td><td>&lt; 0.5</td></tr><tr><td>PDF+&alpha;<sub>s</sub></td><td>Normalization</td><td>0.7</td><td>&lt; 0.5</td></tr><tr><td>Spurious signal</td><td>Normalization</td><td>5.5</td><td>5.4</td></tr></tbody></table>
# 8 Conclusions

Searches for non-resonant and resonant Higgs boson pair production are performed in the $b\bar{b}\gamma\gamma$ final state using 139 fb$^{-1}$ of $pp$ collision data at a center-of-mass energy of 13 TeV. No significant excess with respect to the Standard Model background expectation is observed. A 95% CL upper limit of 130 fb is set on the $pp \to HH$ non-resonant production cross section, where the expected limit is 180 fb. The observed (expected) limit corresponds to 4.1 (5.5) times the cross section predicted by the Standard Model. Constraints on the Higgs boson self-coupling are also derived and limits of $-1.5 < \kappa_\lambda < 6.7$ are obtained, where $-2.4 < \kappa_\lambda < 7.7$ are expected. The constraints on $\kappa_\lambda$ are obtained over an expected hypothesis excluding $pp \to HH$ production. For the resonant production of a scalar particle $X \to HH \to b\bar{b}\gamma\gamma$ upper limits on the production cross section are obtained for the narrow width hypothesis as a function of $m_X$. The observed (expected) upper limits are in the range 610–47 fb (360–43 fb) for 251 GeV $\le m_X \le 1000$ GeV. Compared to the previous ATLAS result based on 36 fb$^{-1}$ of 13 TeV $pp$ collisions, the present analysis extends the data set by more than a factor of four, incorporates a categorization based on $m^*_{b\bar{b}\gamma\gamma}$ and multivariate event selections, and expands the analyzed mass range of the resonance search to lower values. The results improve the previous ATLAS limits on the $HH \to b\bar{b}\gamma\gamma$ production cross section by up to a factor of five and the allowed $\kappa_\lambda$ range shrinks by a factor of $\sim 2$. For the resonant search, the expected limit on the cross section improves by a factor of two to three depending on the $m_X$ value. Of those improvements, a factor of two arises from the increase of integrated luminosity, while the additional improvement can be attributed to the employment of multivariate techniques, more precise object reconstruction and calibration and, for the non-resonant search, to the categorization based on $m^*_{b\bar{b}\gamma\gamma}$. The present analysis improves as well as the constraints set by the ATLAS combination of searches for $HH$ production with up to 36 fb$^{-1}$ of 13 TeV data.
<figure><img src="image_16.png" /><figcaption>Figure 15: Candidate $HH \to b\bar{b}γγ$ event of the high mass BDT tight category. Tracks (in green) with $p_T > 2$ GeV are shown along with two b-jets (red cones) with $p_T$ of 153 GeV and 81 GeV and a $m_{b\bar{b}}$ invariant mass of 113 GeV, two photons (cyan towers) with a transverse energy $E_T$ of 144 GeV and 96 GeV and a $m_{γγ}$ invariant mass of 123 GeV. The $m_{*_{b\bar{b}γγ}}$ invariant mass is 625 GeV.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16: Candidate $HH \to b\bar{b}yy$ event of the high mass BDT tight category. Tracks (in green) with $p_T > 2$ GeV are shown along with two b-jets (red cones) with $p_T$ of 153 GeV and 81 GeV and a $m_{b\bar{b}}$ invariant mass of 113 GeV, two photons (cyan towers) with a transverse energy $E_T$ of 144 GeV and 96 GeV and a $m_{yy}$ invariant mass of 123 GeV. The $m_{b\bar{b}yy}^*$ invariant mass is 625 GeV.</figcaption></figure>
References

[1] ATLAS Collaboration, *Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC*, Phys. Lett. B **716** (2012) 1, arXiv: 1207.7214 [hep-ex] (cit. on p. 2).

[2] CMS Collaboration, *Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC*, Phys. Lett. B **716** (2012) 30, arXiv: 1207.7235 [hep-ex] (cit. on p. 2).

[3] F. Englert and R. Brout, *Broken Symmetry and the Mass of Gauge Vector Mesons*, Phys. Rev. Lett. **13** (1964) 321, ed. by J. C. Taylor (cit. on p. 2).

[4] P. W. Higgs, *Broken symmetries, massless particles and gauge fields*, Phys. Lett. **12** (1964) 132 (cit. on p. 2).

[5] P. W. Higgs, *Broken Symmetries and the Masses of Gauge Bosons*, Phys. Rev. Lett. **13** (1964) 508, ed. by J. C. Taylor (cit. on p. 2).

[6] G. S. Guralnik, C. R. Hagen, and T. W. B. Kibble, *Global Conservation Laws and Massless Particles*, Phys. Rev. Lett. **13** (1964) 585, ed. by J. C. Taylor (cit. on p. 2).

[7] P. W. Higgs, *Spontaneous Symmetry Breakdown without Massless Bosons*, Phys. Rev. **145** (1966) 1156 (cit. on p. 2).

[8] T. W. B. Kibble, *Symmetry breaking in nonAbelian gauge theories*, Phys. Rev. **155** (1967) 1554, ed. by J. C. Taylor (cit. on p. 2).

[9] P. Huang, A. J. Long, and L.-T. Wang, *Probing the Electroweak Phase Transition with Higgs Factories and Gravitational Waves*, Phys. Rev. D **94** (2016) 075008, arXiv: 1608.06619 [hep-ph] (cit. on p. 2).

[10] S. Dawson, S. Dittmaier, and M. Spira, *Neutral Higgs-boson pair production at hadron colliders: QCD corrections*, Phys. Rev. D **58** (1998) 115012, arXiv: hep-ph/9805244 (cit. on p. 2).

[11] S. Borowka et al., *Higgs Boson Pair Production in Gluon Fusion at Next-to-Leading Order with Full Top-Quark Mass Dependence*, Phys. Rev. Lett. **117** (2016) 012001, [Erratum: Phys.Rev.Lett. 117, 079901 (2016)], arXiv: 1604.06447 [hep-ph] (cit. on pp. 2, 5), Erratum: Phys.Rev.Lett. **117** (2016) 079901.

[12] D. de Florian and J. Mazzitelli, *Higgs Boson Pair Production at Next-to-Next-to-Leading Order in QCD*, Phys. Rev. Lett. **111** (2013) 201801, arXiv: 1309.6594 [hep-ph] (cit. on p. 2).

[13] J. Alison et al., *Higgs boson potential at colliders: Status and perspectives*, Rev. Phys. **5** (2020) 100045, ed. by B. Di Micco, M. Gouzevitch, J. Mazzitelli, and C. Vernieri, arXiv: 1910.00012 [hep-ph] (cit. on pp. 2, 18).

[14] D. Y. Shao, C. S. Li, H. T. Li, and J. Wang, *Threshold resummation effects in Higgs boson pair production at the LHC*, JHEP **07** (2013) 169, arXiv: 1301.1245 [hep-ph] (cit. on p. 2).

[15] D. de Florian and J. Mazzitelli, *Higgs pair production at next-to-next-to-leading logarithmic accuracy at the LHC*, JHEP **09** (2015) 053, arXiv: 1505.07122 [hep-ph] (cit. on p. 2).
[16] J. Baglio et al., Gluon fusion into Higgs pairs at NLO QCD and the top mass scheme, Eur. Phys. J. C 79 (2019) 459, arXiv: 1811.05692 [hep-ph] (cit. on p. 2).

[17] M. Grazzini et al., Higgs boson pair production at NNLO with top quark mass effects, JHEP 05 (2018) 059, arXiv: 1803.02463 [hep-ph] (cit. on p. 2).

[18] R. Frederix et al., Higgs pair production at the LHC with NLO and parton-shower effects, Phys. Lett. B 732 (2014) 142, arXiv: 1401.7340 [hep-ph] (cit. on p. 2).

[19] G. D. Kribs and A. Martin, Enhanced di-Higgs Production through Light Colored Scalars, Phys. Rev. D 86 (2012) 095023, arXiv: 1207.4496 [hep-ph] (cit. on p. 2).

[20] L. Di Luzio, R. Gröber, and M. Spannowsky, Maxi-sizing the trilinear Higgs self-coupling: how large could it be? Eur. Phys. J. C 77 (2017) 788, arXiv: 1704.02311 [hep-ph] (cit. on p. 2).

[21] G. D. Kribs, A. Maier, H. Rzehak, M. Spannowsky, and P. Waite, Electroweak oblique parameters as a probe of the trilinear Higgs boson self-interaction, Phys. Rev. D 95 (2017) 093004, arXiv: 1702.07678 [hep-ph] (cit. on p. 2).

[22] T. D. Lee, A Theory of Spontaneous T Violation, Phys. Rev. D 8 (1973) 1226, ed. by G. Feinberg (cit. on p. 3).

[23] S. Dimopoulos and H. Georgi, Softly Broken Supersymmetry and SU(5), Nucl. Phys. B 193 (1981) 150 (cit. on p. 3).

[24] Z. Chacko, Y. Nomura, M. Papucci, and G. Perez, Natural Little Hierarchy from a Partially Goldstone Twin Higgs, JHEP 01 (2006) 126, arXiv: hep-ph/0510273 (cit. on p. 3).

[25] J. Mrazek et al., The Other Natural Two Higgs Doublet Model, Nucl. Phys. B 853 (2011) 1, arXiv: 1105.5403 [hep-ph] (cit. on p. 3).

[26] L. Randall and R. Sundrum, Large Mass Hierarchy from a Small Extra Dimension, Phys. Rev. Lett. 83 (1999) 3370, arXiv: hep-ph/9905221 (cit. on p. 3).

[27] ATLAS Collaboration, Search for Higgs boson pair production in the $\gamma\gamma b\bar{b}$ final state with 13 TeV pp collision data collected by the ATLAS experiment, JHEP 11 (2018) 040, arXiv: 1807.04873 [hep-ex] (cit. on p. 3).

[28] ATLAS Collaboration, Combination of searches for Higgs boson pairs in pp collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector, Phys. Lett. B 800 (2020) 135103, arXiv: 1906.02025 [hep-ex] (cit. on pp. 3, 4).

[29] CMS Collaboration, Search for nonresonant Higgs boson pair production in final states with two bottom quarks and two photons in proton-proton collisions at $\sqrt{s} = 13$ TeV, (2020), arXiv: 2011.12373 [hep-ex] (cit. on p. 4).

[30] CMS Collaboration, Search for Higgs boson pair production in the $\gamma\gamma b\bar{b}$ final state in pp collisions at $\sqrt{s} = 13$ TeV, Phys. Lett. B 788 (2019) 7, arXiv: 1806.00408 [hep-ex] (cit. on p. 4).

[31] ATLAS Collaboration, The ATLAS Experiment at the CERN Large Hadron Collider, JINST 3 (2008) S08003 (cit. on p. 4).

[32] ATLAS Collaboration, ATLAS Insertable B-Layer: Technical Design Report, ATLAS-TDR-19; CERN-LHCC-2010-013, 2010, URL: https://cds.cern.ch/record/1291633 (cit. on p. 4).
[33] B. Abbott et al., Production and integration of the ATLAS Insertable B-Layer, JINST 13 (2018) T05008, arXiv: 1803.00844 [physics.ins-det] (cit. on p. 4).

[34] ATLAS Collaboration, Performance of the ATLAS trigger system in 2015, Eur. Phys. J. C 77 (2017) 317, arXiv: 1611.09661 [hep-ex] (cit. on p. 4).

[35] ATLAS Collaboration, ATLAS data quality operations and performance for 2015–2018 data-taking, JINST 15 (2020) P04003, arXiv: 1911.04632 [physics.ins-det] (cit. on p. 5).

[36] G. Avoni et al., The new LUCID-2 detector for luminosity measurement and monitoring in ATLAS, JINST 13 (2018) P07017 (cit. on pp. 5, 17).

[37] ATLAS Collaboration, Luminosity determination in pp collisions at $\sqrt{s}$ = 13 TeV using the ATLAS detector at the LHC, ATLAS-CONF-2019-021, 2019, URL: https://cds.cern.ch/record/2677054 (cit. on pp. 5, 17).

[38] ATLAS Collaboration, Emulating the impact of additional proton-proton interactions in the ATLAS simulation by pre-sampling sets of inelastic Monte Carlo events, (2021), arXiv: 2102.09495 [hep-ex] (cit. on p. 5).

[39] P. Nason and C. Oleari, NLO Higgs boson production via vector-boson fusion matched with shower in POWHEG, JHEP 02 (2010) 037, arXiv: 0911.5299 [hep-ph] (cit. on pp. 5, 6).

[40] G. Heinrich, S. P. Jones, M. Kerner, G. Luisoni, and E. Vryonidou, NLO predictions for Higgs boson pair production with full top quark mass dependence matched to parton showers, JHEP 08 (2017) 088, arXiv: 1703.09252 [hep-ph] (cit. on p. 5).

[41] G. Heinrich, S. Jones, M. Kerner, G. Luisoni, and L. Scyboz, Probing the trilinear Higgs boson coupling in di-Higgs production at NLO QCD including parton shower effects, JHEP 06 (2019) 066, arXiv: 1903.08137 [hep-ph] (cit. on p. 5).

[42] J. Butterworth et al., PDF4LHC recommendations for LHC Run II, J. Phys. G 43 (2016) 023001, arXiv: 1510.03865 [hep-ph] (cit. on pp. 5, 6).

[43] S. Amoroso et al., “Les Houches 2019: Physics at TeV Colliders: Standard Model Working Group Report,” 11th Les Houches Workshop on Physics at TeV Colliders: PhysTeV Les Houches, 2020, arXiv: 2003.01700 [hep-ph] (cit. on p. 5).

[44] J. Alwall et al., The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP 07 (2014) 079, arXiv: 1405.0301 [hep-ph] (cit. on pp. 5, 6).

[45] F. Bishara, R. Contino, and J. Rojo, Higgs pair production in vector-boson fusion at the LHC and beyond, Eur. Phys. J. C 77 (2017) 481, arXiv: 1611.03860 [hep-ph] (cit. on p. 5).

[46] R. D. Ball et al., Parton distributions with LHC data, Nucl. Phys. B867 (2013) 244, arXiv: 1207.1303 [hep-ph] (cit. on p. 5).

[47] L.-S. Ling et al., NNLO QCD corrections to Higgs pair production via vector boson fusion at hadron colliders, Phys. Rev. D 89 (2014) 073001, arXiv: 1401.7754 [hep-ph] (cit. on p. 5).
[48] F. A. Dreyer and A. Karlberg, *Vector-Boson Fusion Higgs Pair Production at N³LO*, Phys. Rev. D **98** (2018) 114016, arXiv: 1811.07906 [hep-ph] (cit. on p. 5).

[49] F. A. Dreyer and A. Karlberg, *Fully differential Vector-Boson Fusion Higgs Pair Production at Next-to-Next-to-Leading Order*, Phys. Rev. D **99** (2019) 074028, arXiv: 1811.07918 [hep-ph] (cit. on p. 5).

[50] M. Bähr et al., *Herwig++ physics and manual*, Eur. Phys. J. C **58** (2008) 639, arXiv: 0803.0883 [hep-ph] (cit. on p. 5).

[51] J. Bellm et al., *Herwig 7.0/Herwig++ 3.0 release note*, Eur. Phys. J. C **76** (2016) 196, arXiv: 1512.01178 [hep-ph] (cit. on p. 5).

[52] ATLAS Collaboration, *Combined measurements of Higgs boson production and decay using up to 80 fb⁻¹ of proton-proton collision data at √s = 13 TeV collected with the ATLAS experiment*, Phys. Rev. D **101** (2020) 012002, arXiv: 1909.02845 [hep-ex] (cit. on pp. 5, 19).

[53] ATLAS and CMS Collaborations, *Combined Measurement of the Higgs Boson Mass in pp Collisions at √s = 7 and 8 TeV with the ATLAS and CMS Experiments*, Phys. Rev. Lett. **114** (2015) 191803, arXiv: 1503.07589 [hep-ex] (cit. on pp. 6, 18).

[54] D. de Florian et al., *Handbook of LHC Higgs Cross Sections: 4. Deciphering the Nature of the Higgs Sector*, (2016), arXiv: 1610.07922 [hep-ph] (cit. on pp. 6, 18).

[55] A. Djouadi, J. Kalinowski, and M. Spira, *HDECAY: A Program for Higgs boson decays in the standard model and its supersymmetric extension*, Comput. Phys. Commun. **108** (1998) 56, arXiv: hep-ph/9704448 (cit. on p. 6).

[56] E. Bothmann et al., *Event generation with Sherpa 2.2*, SciPost Phys. **7** (2019) 034, arXiv: 1905.09127 [hep-ph] (cit. on p. 6).

[57] T. Gleisberg and S. Höhe, *Comix, a new matrix element generator*, JHEP **12** (2008) 039, arXiv: 0808.3674 [hep-ph] (cit. on p. 6).

[58] F. Buccioni et al., *OpenLoops 2*, Eur. Phys. J. C **79** (2019) 866, arXiv: 1907.13071 [hep-ph] (cit. on p. 6).

[59] F. Cascioli, P. Maierhöfer, and S. Pozzorini, *Scattering Amplitudes with Open Loops*, Phys. Rev. Lett. **108** (2012) 111601, arXiv: 1111.5206 [hep-ph] (cit. on p. 6).

[60] A. Denner, S. Dittmaier, and L. Hofer, *COLLIER: A fortran-based complex one-loop library in extended regularizations*, Comput. Phys. Commun. **212** (2017) 220, arXiv: 1604.06792 [hep-ph] (cit. on p. 6).

[61] S. Schumann and F. Krauss, *A parton shower algorithm based on Catani–Seymour dipole factorisation*, JHEP **03** (2008) 038, arXiv: 0709.1027 [hep-ph] (cit. on p. 6).

[62] S. Höhe, F. Krauss, M. Schönherr, and F. Siegert, *A critical appraisal of NLO+PS matching methods*, JHEP **09** (2012) 049, arXiv: 1111.1220 [hep-ph] (cit. on p. 6).

[63] S. Höhe, F. Krauss, M. Schönherr, and F. Siegert, *QCD matrix elements + parton showers. The NLO case*, JHEP **04** (2013) 027, arXiv: 1207.5030 [hep-ph] (cit. on p. 6).
[64] F. Siegert, *A practical guide to event generation for prompt photon production with Sherpa*, J. Phys. G 44 (2017) 044007, arXiv: 1611.07226 [hep-ph] (cit. on p. 6).

[65] K. Hamilton, P. Nason, and G. Zanderighi, *MINLO: Multi-Scale Improved NLO*, JHEP 10 (2012) 155, arXiv: 1206.3572 [hep-ph] (cit. on p. 6).

[66] J. M. Campbell et al., *NLO Higgs Boson Production Plus One and Two Jets Using the POWHEG BOX*, MadGraph4 and MCFM, JHEP 07 (2012) 092, arXiv: 1202.5475 [hep-ph] (cit. on p. 6).

[67] K. Hamilton, P. Nason, C. Oleari, and G. Zanderighi, *Merging H/W/Z + 0 and 1 jet at NLO with no merging scale: a path to parton shower + NNLO matching*, JHEP 05 (2013) 082, arXiv: 1212.4504 [hep-ph] (cit. on p. 6).

[68] G. Bozzi, S. Catani, D. de Florian, and M. Grazzini, *Transverse-momentum resummation and the spectrum of the Higgs boson at the LHC*, Nucl. Phys. B 737 (2006) 73, arXiv: hep-ph/0508068 [hep-ph] (cit. on p. 6).

[69] D. de Florian, G. Ferrera, M. Grazzini, and D. Tommasini, *Transverse-momentum resummation: Higgs boson production at the Tevatron and the LHC*, JHEP 11 (2011) 064, arXiv: 1109.2109 [hep-ph] (cit. on p. 6).

[70] T. Sjöstrand et al., *An introduction to PYTHIA 8.2*, Comput. Phys. Commun. 191 (2015) 159, arXiv: 1410.3012 [hep-ph] (cit. on p. 6).

[71] ATLAS Collaboration, *Measurement of the Z/γ*\* boson transverse momentum distribution in pp collisions at √s = 7 TeV with the ATLAS detector*, JHEP 09 (2014) 145, arXiv: 1406.3660 [hep-ex] (cit. on p. 6).

[72] K. Hamilton, P. Nason, E. Re, and G. Zanderighi, *NNLOPS simulation of Higgs boson production*, JHEP 10 (2013) 222, arXiv: 1309.0017 [hep-ph] (cit. on p. 6).

[73] P. Nason, *A new method for combining NLO QCD with shower Monte Carlo algorithms*, JHEP 11 (2004) 040, arXiv: hep-ph/0409146 (cit. on p. 6).

[74] S. Frixione, P. Nason, and C. Oleari, *Matching NLO QCD computations with parton shower simulations: the POWHEG method*, JHEP 11 (2007) 070, arXiv: 0709.2092 [hep-ph] (cit. on p. 6).

[75] S. Alioli, P. Nason, C. Oleari, and E. Re, *A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX*, JHEP 06 (2010) 043, arXiv: 1002.2581 [hep-ph] (cit. on p. 6).

[76] K. Mimasu, V. Sanz, and C. Williams, *Higher order QCD predictions for associated Higgs production with anomalous couplings to gauge bosons*, JHEP 08 (2016) 039, arXiv: 1512.02572 [hep-ph] (cit. on p. 6).

[77] G. Luisoni, P. Nason, C. Oleari, and F. Tramontano, $HW^+/HZ + 0$ and 1 jet at NLO with the POWHEG BOX interfaced to GoSam and their merging within MiNLO, JHEP 10 (2013) 083, arXiv: 1306.2542 [hep-ph] (cit. on p. 6).

[78] H. B. Hartanto, B. Jäger, L. Reina, and D. Wackeroth, *Higgs boson production in association with top quarks in the POWHEG BOX*, Phys. Rev. D 91 (2015) 094003, arXiv: 1501.04498 [hep-ph] (cit. on p. 6).

[79] S. Frixione, P. Nason, and G. Ridolfi, *A positive-weight next-to-leading-order Monte Carlo for heavy flavour hadroproduction*, JHEP 09 (2007) 126, arXiv: 0707.3088 [hep-ph] (cit. on p. 6).
[80] R. D. Ball et al., Parton distributions for the LHC run II, JHEP 04 (2015) 040, arXiv: 1410.8849 [hep-ph] (cit. on p. 6).

[81] ATLAS Collaboration, ATLAS Pythia 8 tunes to 7 TeV data, ATL-PHYS-PUB-2014-021, 2014, URL: https://cds.cern.ch/record/1966419 (cit. on p. 6).

[82] ATLAS Collaboration, The Pythia 8 A3 tune description of ATLAS minimum bias and inelastic measurements incorporating the Donnachie-Landshoff diffractive model, ATL-PHYS-PUB-2016-017, 2016, URL: https://cds.cern.ch/record/2206965 (cit. on p. 6).

[83] ATLAS Collaboration, The ATLAS Simulation Infrastructure, Eur. Phys. J. C 70 (2010) 823, arXiv: 1005.4568 [physics.ins-det] (cit. on p. 6).

[84] GEANT4 Collaboration, S. Agostinelli, et al., GEANT4 – a simulation toolkit, Nucl. Instrum. Meth. A 506 (2003) 250 (cit. on p. 6).

[85] ATLAS Collaboration, The simulation principle and performance of the ATLAS fast calorimeter simulation FastCaloSim, (2010), ATL-PHYS-PUB-2010-013, URL: http://cdsweb.cern.ch/record/1300517 (cit. on p. 6).

[86] ATLAS Collaboration, Electron and photon reconstruction and performance in ATLAS using a dynamical, topological cell clustering-based approach, ATL-PHYS-PUB-2017-022, 2017, URL: https://cds.cern.ch/record/2298955 (cit. on p. 7).

[87] ATLAS Collaboration, Electron and photon performance measurements with the ATLAS detector using the 2015–2017 LHC proton–proton collision data, JINST 14 (2019) P12006, arXiv: 1908.00005 [hep-ex] (cit. on pp. 7, 9, 17).

[88] ATLAS Collaboration, Measurement of Higgs boson production in the diphoton decay channel in pp collisions at center-of-mass energies of 7 and 8 TeV with the ATLAS detector, Phys. Rev. D 90 (2014) 112015, arXiv: 1408.7084 [hep-ex] (cit. on pp. 7, 16, 17).

[89] ATLAS Collaboration, Muon reconstruction performance of the ATLAS detector in proton–proton collision data at $\sqrt{s} = 13$ TeV, Eur. Phys. J. C 76 (2016) 292, arXiv: 1603.05598 [hep-ex] (cit. on p. 7).

[90] ATLAS Collaboration, Muon reconstruction and identification efficiency in ATLAS using the full Run 2 pp collision data set at $\sqrt{s} = 13$ TeV, (2020), arXiv: 2012.00578 [hep-ex] (cit. on p. 7).

[91] ATLAS Collaboration, Jet reconstruction and performance using particle flow with the ATLAS Detector, Eur. Phys. J. C 77 (2017) 466, arXiv: 1703.10485 [hep-ex] (cit. on p. 8).

[92] M. Cacciari, G. P. Salam, and G. Soyez, FastJet User Manual, Eur. Phys. J. C 72 (2012) 1896, arXiv: 1111.6097 [hep-ph] (cit. on p. 8).

[93] M. Cacciari, G. P. Salam, and G. Soyez, The anti-$k_t$ jet clustering algorithm, JHEP 04 (2008) 063, arXiv: 0802.1189 [hep-ph] (cit. on p. 8).

[94] ATLAS Collaboration, Tagging and suppression of pileup jets with the ATLAS detector, ATLAS-CONF-2014-018, 2014, URL: https://cds.cern.ch/record/1700870 (cit. on p. 8).

[95] ATLAS Collaboration, ATLAS b-jet identification performance and efficiency measurement with $t\bar{t}$ events in pp collisions at $\sqrt{s} = 13$ TeV, Eur. Phys. J. C 79 (2019) 970, arXiv: 1907.05120 [hep-ex] (cit. on pp. 8, 18).
[96] ATLAS Collaboration, *Identification of Jets Containing b-Hadrons with Recurrent Neural Networks at the ATLAS Experiment*, ATL-PHYS-PUB-2017-003, 2017, URL: https://cds.cern.ch/record/2255226 (cit. on p. 8).

[97] ATLAS Collaboration, *Optimisation and performance studies of the ATLAS b-tagging algorithms for the 2017-18 LHC run*, ATL-PHYS-PUB-2017-013, 2017, URL: https://cds.cern.ch/record/2273281 (cit. on p. 8).

[98] ATLAS Collaboration, *Performance of b-jet identification in the ATLAS experiment*, JINST 11 (2016) P04008, arXiv: 1512.01094 [hep-ex] (cit. on p. 8).

[99] ATLAS Collaboration, *Evidence for the $H \to b\bar{b}$ decay with the ATLAS detector*, JHEP 12 (2017) 024, arXiv: 1708.03299 [hep-ex] (cit. on p. 8).

[100] ATLAS Collaboration, *Performance of missing transverse momentum reconstruction with the ATLAS detector using proton-proton collisions at $\sqrt{s}$ = 13 TeV*, Eur. Phys. J. C 78 (2018) 903, arXiv: 1802.08168 [hep-ex] (cit. on pp. 8, 18).

[101] ATLAS Collaboration, *Performance of electron and photon triggers in ATLAS during LHC Run 2*, Eur. Phys. J. C 80 (2020) 47, arXiv: 1909.00761 [hep-ex] (cit. on pp. 9, 17).

[102] ATLAS Collaboration, *Search for pair production of Higgs bosons in the $b\bar{b}b$ final state using proton-proton collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector*, JHEP 01 (2019) 030, arXiv: 1804.06174 [hep-ex] (cit. on p. 9).

[103] T. Chen and C. Guestrin, “XGBoost: A Scalable Tree Boosting System,” KDD 16, (2016), 785, arXiv: 1603.02754 [cs.LG] (cit. on p. 9).

[104] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, *Asymptotic formulae for likelihood-based tests of new physics*, Eur. Phys. J. C 71 (2011) 1554, arXiv: 1007.1727 [physics.data-an] (cit. on pp. 11, 19), Erratum: Eur. Phys. J. C 73 (2013) 2501.

[105] A. Hoecker et al., *TMVA - Toolkit for Multivariate Data Analysis*, (2007), arXiv: physics/0703039 (cit. on p. 12).

[106] ATLAS Collaboration, *Measurement of isolated-photon pair production in pp collisions at $\sqrt{s}$ = 7 TeV with the ATLAS detector*, JHEP 01 (2013) 086, arXiv: 1211.1913 [hep-ex] (cit. on p. 13).

[107] M. Oreglia, *A Study of the Reactions $\psi' \to \gamma\gamma\psi$, 1980*, URL: https://www.slac.stanford.edu/cgi-wrap/getdoc/slac-r-236.pdf (cit. on p. 16).

[108] ATLAS Collaboration, *Recommendations for the Modeling of Smooth Backgrounds*, ATL-PHYS-PUB-2020-028, 2020, URL: https://cds.cern.ch/record/2743717 (cit. on p. 17).

[109] A. Wald, *Sequential Tests of Statistical Hypotheses*, Ann. Math. Statist. 16 (1945) 117 (cit. on p. 17).

[110] ATLAS Collaboration, *Search for resonances in diphoton events at $\sqrt{s}$ = 13 TeV with the ATLAS detector*, JHEP 09 (2016) 001, arXiv: 1606.03833 [hep-ex] (cit. on p. 17).

[111] ATLAS Collaboration, *Measurement of b-tagging efficiency of c-jets in $t\bar{t}$ events using a likelihood approach with the ATLAS detector*, ATLAS-CONF-2018-001, 2018, URL: https://cds.cern.ch/record/2306649 (cit. on p. 18).
[112] ATLAS Collaboration, *Calibration of light-flavour b-jet mistagging rates using ATLAS proton-proton collision data at $\sqrt{s}$ = 13 TeV*, ATLAS-CONF-2018-006, 2018, URL: https://cds.cern.ch/record/2314418 (cit. on p. 18).

[113] ATLAS Collaboration, *Measurements of inclusive and differential fiducial cross-sections of $t\bar{t}$ production with additional heavy-flavour jets in proton-proton collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector*, JHEP 04 (2019) 046, arXiv: 1811.12113 [hep-ex] (cit. on p. 18).

[114] ATLAS Collaboration, *Study of heavy-flavor quarks produced in association with top quark pairs at $\sqrt{s}$ = 7 TeV using the ATLAS detector*, Phys. Rev. D 89 (2014) 072012, arXiv: 1304.6386 [hep-ex] (cit. on p. 18).

[115] ATLAS Collaboration, *Measurement of the cross-section for W boson production in association with b-jets in pp collisions at $\sqrt{s}$ = 7 TeV with the ATLAS detector*, JHEP 06 (2013) 084, arXiv: 1302.2929 [hep-ex] (cit. on p. 18).

[116] A. L. Read, *Presentation of search results: the CLS technique*, J. Phys. G 28 (2002) 2693 (cit. on p. 19).