Search for Higgs boson pair production in the two bottom quarks plus two photons final state in pp collisions at

Searches are performed for nonresonant and resonant di-Higgs boson production in the b ¯ b γγ final state. The dataset 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 above the expected background is found and upper limits on the di-Higgs boson production cross sections are set. A 95% confidence-level upper limit of 4.2 times the cross section predicted by the Standard Model is set on pp → HH nonresonant production, where the expected limit is 5.7 times the Standard Model predicted value. The expected constraints are obtained for a background hypothesis excluding pp → HH production. The observed (expected) constraints on the Higgs boson trilinear coupling modifier κ λ are determined to be ½ − 1 . 5 ; 6 . 7 (cid:2) ð½ − 2 . 4 ; 7 . 7 (cid:2)Þ at 95% confidence level, where the expected constraints on κ λ are obtained excluding pp → HH production from the background hypothesis. For resonant production of a new hypothetical scalar particle X ( X → HH → b ¯ b γγ ), limits on the cross section for pp → X → HH are presented in the narrow-width approximation as a function of m X in the range 251 GeV ≤ m X ≤ 1000 GeV. The observed (expected) limits on the cross section for pp → X → HH range from 640 fb to 44 fb (391 fb to 46 fb) over the considered mass range.


Search for Higgs boson pair production in the two bottom quarks plus two photons final state in collisions at √ = 13 TeV with the ATLAS detector
The ATLAS Collaboration Searches are performed for nonresonant and resonant di-Higgs boson production in thē 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 above the expected background is found and upper limits on the di-Higgs boson production cross sections are set. A 95% confidence-level upper limit of 4.2 times the cross section predicted by the Standard Model is set on → nonresonant production, where the expected limit is 5.7 times the Standard Model predicted value. The expected constraints are obtained for a background hypothesis excluding → production. The observed (expected) constraints on the Higgs boson trilinear coupling modifier are determined to be [−1.5, 6.7] ( [−2.4, 7.7]) at 95% confidence level, where the expected constraints on are obtained excluding → production from the background hypothesis. For resonant production of a new hypothetical scalar particle ( → →¯), limits on the cross section for → → are presented in the narrow-width approximation as a function of in the range 251 ≤ ≤ 1000 GeV. The observed (expected) limits on the cross section for → → range from 640 fb to 44 fb (391 fb to 46 fb) over the considered mass range.

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][4][5][6][7][8]. The Higgs boson self-coupling provides information about the structure of the Higgs potential. It is possible to directly probe the self-coupling of the Higgs boson by studying Higgs boson pair ( ) production. Furthermore, any deviation of the Higgs boson pair production rate from the Standard Model (SM) prediction would point to new physics beyond the Standard Model (BSM) and may be within the sensitivity reach of the proton-proton ( ) collision data collected at √ = 13 TeV during 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][11][12]. For 13 TeV collisions and a Higgs boson mass = 125.09 GeV [13], the ggF cross section, calculated at next-to-next-to-leading-order (NNLO) accuracy in the finite top-quark mass approximation (FTapprox), is (ggF) = 31.02 +2.2% −5.0% (Scale) +4% −18% ( top ) ± 3.0% ( s +PDF) fb [14][15][16][17], where "Scale" represents the uncertainty due to the finite order of the quantum chromodynamics (QCD) calculation, " top " the uncertainty related to the top-quark mass scheme [17,18] which is added linearly to the Scale uncertainty, and " s +PDF" the effect of uncertainties in the strong coupling constant and parton distribution functions. The di-Higgs vector-boson fusion (VBF) production cross section, calculated at next-to-next-to-next-toleading order (N3LO) for = 125.09 GeV, is (VBF) = 1.72 +0.03% −0.04% (Scale) ±2.1% ( s +PDF) fb [14], which is one order of magnitude lower than the cross section of the ggF process. The VBF production mode provides the analysis with additional sensitivity to the Higgs trilinear coupling , as shown in Figure 2. Both the ggF and VBF production modes of Higgs boson pairs are considered as signal modes in this paper. The other production modes have lower cross sections [19] and are neglected. Nonresonant enhancements to the Higgs boson pair cross section can originate either from loop corrections involving new particles, such as light, colored scalars [20], or from non-SM couplings between the Higgs boson and other SM particles. The nonresonant production cross section can also be altered by the trilinear self-coupling, , being different from the SM prediction, as discussed in Refs. [21,22]. Such an effect can be captured by a scale factor defined as = / SM , where SM is the SM value of the parameter.
In addition to the nonresonant 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 existence of heavy scalar particles that can decay into a pair of Higgs bosons. Such theories include models with two Higgs doublets [23], such as the minimal supersymmetric extension of the SM [24], twin Higgs models [25] and composite Higgs models [26], adding a second complex scalar doublet to the Higgs sector. Alternatively, the Randall-Sundrum model of warped extra dimensions [27] predicts spin-0 radions that could couple to a Higgs boson pair. This paper presents a search for di-Higgs production in the¯final state, including dedicated assessments of nonresonant and resonant contributions. The analysis considers the full Run 2 data set of 139 fb −1 at 13 TeV. For both the nonresonant 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 diphoton invariant mass, . For the nonresonant search, data are divided into different categories based on the four-body invariant mass to target different ranges. The resonant search focuses on probing the existence of a narrow-width scalar particle in the mass range 251 < < 1000 GeV decaying into a pair of Higgs bosons. The selection criteria depend on the mass of the probed scalar particle. The main background processes are diphoton-plus-jets production and processes where a Higgs boson is produced and decays into a pair of photons. In the context of the resonant search, nonresonant production is considered as a background.
Previous results from 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 they were found to be consistent with SM expectations within uncertainties [28]. The search for nonresonant enhancements of Higgs boson pair production set an observed (expected) 95% confidence level (CL) upper limit on the cross section of 0.73 (0.93) pb, corresponding to 22 (28) times the SM prediction. The Higgs trilinear coupling was constrained to the range −8.2 < < 13.2 at 95% CL (−8.3 < < 13.2 expected). A previous combination of searches for pair production performed by the ATLAS Collaboration with up to 36 fb −1 of 13 TeV data provided constraints of −5.0 < < 12.0 at 95% CL [29]. A search for enhancements due to the decay of a narrow-width scalar particle was performed in this channel as well, and results were given as a function of the resonance mass, . The observed (expected) limits were between 1.1 pb and 0.12 pb over the range 260 < < 1000 GeV. A combination of searches by the ATLAS Collaboration with up to 36 fb −1 of 13 TeV data for a narrow-width scalar resonance decaying into a pair was performed and provided upper limits between 851 fb and 4.6 fb over the range 260 < < 3000 GeV [29]. The CMS Collaboration also set observed (expected) upper limits of 7.7 (5.2) times the SM prediction at 95% CL on nonresonant enhancements of Higgs boson pair production in the¯final state with 137 fb −1 of 13 TeV data [30]. In the same final state, 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 < < 900 GeV with 35.9 fb −1 of 13 TeV data [31].

The ATLAS detector
The ATLAS detector [32] at the LHC covers nearly the entire solid angle around the collision point. 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 | | < 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 installed before Run 2 [33,34]. It is followed by the silicon microstrip tracker, 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 | | = 2.0. The TRT also provides electron identification information.
The calorimeter system covers the pseudorapidity range | | < 4.9. Within the region | | < 3.2, electromagnetic calorimetry is provided by barrel and endcap high-granularity lead/liquid-argon (LAr) calorimeters, with an additional thin LAr presampler covering | | < 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 | | < 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 | | < 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 | | < 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 [35]. The first-level trigger accepts events from the 40 MHz bunch crossings at a rate below 100 kHz, which the high-level trigger reduces in order to record events at rate of about 1 kHz. An extensive software suite [36] is used in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.

Data and simulation samples
This analysis uses collision data collected by the ATLAS experiment from 2015 to 2018 with proton beams colliding at a center-of-mass energy of √ = 13 TeV. After data quality requirements [37] the full 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 -axis along the beam pipe. The -axis points from the IP to the center of the LHC ring, and the -axis points upwards. Cylindrical coordinates ( , ) are used in the transverse plane, being the azimuthal angle around the -axis. The pseudorapidity is defined in terms of the polar angle as = − ln tan( /2). The rapidity is defined in terms of the energy, the momentum and the polar angle : = Monte Carlo (MC) simulations are available for the signal as well as most background processes as detailed in the rest of this section. The reducible backgrounds from final states with jets wrongly identified as photons ( -jet and dĳet backgrounds) are, however, estimated using a data-driven technique detailed in Section 4.3.
Events from ggF nonresonant production were generated at next-to-leading-order accuracy in QCD with finite top-quark mass in both the real and virtual corrections (NLO FT) [11], using the P B v2 [41] generator in the finite top-quark mass approximation [42,43] with the PDF4LHC15 parton distribution function (PDF) set [44]. The P 8.244 generator was used for parton showering, hadronization and underlying-event simulation. H 7.1.6 was used as an alternative generator to calculate the theory uncertainty from the parton shower. Samples were generated for coupling modifier values = 1 and 10.
For ggF nonresonant production, a reweighting method based on the di-Higgs invariant mass is used to provide predictions on the cross section at different values, starting from the existing = 1 sample. The reweighting method derives the scale factors as a function of in bins of by performing a linear combination of samples generated at different values [45]. Histograms of the truth distribution are produced for each sample and the distributions of the other relevant kinematic variables are obtained applying an event-per-event weight based on the ratio between the binned distribution for the targeted by the binned distribution for the SM ( = 1) sample.
This method was validated by comparing the event yields and the distributions of the relevant Higgs boson kinematic variables, including the variable, of the sample generated with = 10 to the sample generated with = 1 and reweighted to = 10. Good agreement is obtained in all categories. A systematic uncertainty in the range of 3%-4% is associated with the reweighting process, based on the maximum differences of signal yields observed in this validation. For each value, the inclusive cross section is normalized according to Ref. [46]. A fit with a second order polynomial to the MC prediction is performed in each analysis category, in order to parametrize the event yields as a function of .  [49] was used in the matrix element, interfaced to P 8.244. The cross section of the VBF process is evaluated at N3LO in QCD [50][51][52], as outlined in Section 1. Since samples were generated at LO for four values of the coupling modifier, = 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 cross section. These samples are used to derive a parameterization of the signal yields in the signal region as a function of by fitting a second-order polynomial to the MC predictions in each analysis category, as described in Section 4.2.2.
The production of a heavy spin-0 resonance via ggF and its decay into a pair of Higgs  ( → and → ),¯, ( and ), and was modeled using the same set of MC samples as described in Ref. [55]. For single Higgs boson production, as well as both nonresonant and resonant di-Higgs production, a Higgs boson mass of 125.09 GeV was assumed [13]. The analysis assumes a branching ratio of 0.227% for the Higgs boson decay into two photons and a branching ratio of 58.2% for the Higgs boson decay into two -quarks [56,57]. The inclusive cross sections of these processes are normalized to the most precise available theoretical values [56].
The +jets process was simulated with the S 2.2.4 [58] generator. QCD NLO-accurate matrix elements for up to one parton, and LO-accurate matrix elements for up to three partons, were calculated with the Comix [59] and O L [60][61][62] libraries. These were calculated in the five-flavor scheme including -quarks in the massless approximation and merged with the S parton shower [63] using the MEPS@NLO prescription [64,65] with a dynamic merging cut [66] of 10 GeV. Within the parton shower, -quarks were then treated as being massive. Finally, events from¯processes were produced with M G 5_ MC@NLO in the four-flavor scheme [47]. The simulation samples used in the analysis are listed in Table 1.

Object selection
Photons are reconstructed from topologically connected clusters [89] 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 conversion vertices or tracks that 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 → + − events [90]. The photon direction is reconstructed using the longitudinal (i.e. shower-depth) segmentation of the calorimeters and constrained to be compatible with the luminous region of the proton beam collisions. Additionally, information about the position of the conversion vertex and about the tracks associated with the conversion vertex is considered in the case of converted photons.
Events are required to have at least one reconstructed collision vertex, defined as a vertex associated with at least two tracks with transverse momentum ( T ) larger than 0.5 GeV. The primary vertex is selected from the reconstructed collision vertices using a neural-network algorithm [91] based on the extrapolated photon trajectories and the tracks associated with 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 misidentification of hadronic jets containing large neutral components, primarily 0 particles, which decay into a pair of highly collimated photons. "Tight" identification criteria, which are tuned for converted and unconverted photons separately, are applied [90].
To further improve the rejection of misidentified photons, two isolation variables are defined to quantify the activity around a photon. Calorimeter-based isolation iso T is defined as the sum of the transverse energy of topological clusters within a cone of size Δ = 0.2 around the photon, correcting for the energy of the photon candidate itself as well as for an average expected pileup contribution. Track-based isolation iso T is defined as the scalar sum of the transverse momenta of all tracks with T > 1 GeV that originate from the primary vertex and are within a cone of Δ = 0.2 around the photon. Isolated photons must have iso T < 0.065 · T and iso T < 0.05 · T (the "loose" working point [90]), where T is the transverse energy of the photon. For isolated photons with transverse momenta between 30 GeV and 250 GeV, the identification efficiency for unconverted and converted photons ranges from 84% to 98% [90].
Electrons are reconstructed from energy deposits measured in the electromagnetic calorimeter which are matched to ID tracks [90]. They are required to satisfy | | < 2.47, excluding the calorimeter transition region 1.37 < | | < 1.52, and have a transverse momentum 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 [90]. Muons are reconstructed from high-quality tracks found in the MS [92]. A matching of these tracks to ID tracks is required in the region | | < 2.5. Muons are required to have | | < 2.7 and T > 10 GeV, and to satisfy a "medium" identification criterion [93]. Both the electrons and muons are matched to the primary vertex via requirements on the tracks' longitudinal and transverse impact parameters, | 0 | and | 0 |, respectively. These requirements are | 0 | sin < 0.5 mm (where is the polar angle of the track) for electrons and muons and | 0 |/ 0 < 5 (3) for electrons (muons).
Reconstructed jets are based on particle-flow objects built from noise-suppressed positive-energy topological clusters in the calorimeter and reconstructed tracks [94]. The anti-algorithm [95,96] with a radius parameter of = 0.4 is used. They are required to have rapidity | | < 4.4 and T > 25 GeV. To further suppress jets produced in concurrent interactions, each jet within the tracking acceptance of | | < 2.4, and with T < 60 GeV, is required to satisfy the "tight" jet-vertex tagger [97] 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 -tagging is based on distinctive features of -hadron decays in terms of the impact parameters of the tracks and the displaced vertices reconstructed in the inner detector [98]. The inputs of the DL1r network also include discriminating variables constructed by a recurrent neural network (RNNIP) [99], which exploits the spatial and kinematic correlations between tracks originating from the same -hadron. For high-T jets, this approach is found to give better performance [100] than a previously used multivariate technique [101]. Operating points are defined by a single selection value on the discriminant output distribution and are chosen to provide a specific -jet efficiency for an inclusive¯MC sample. The -tagging requirements result in an efficiency of 77% for jets containing -hadrons ( -tagged jets), and the misidentification rate is 1/130 (1/4.9) for light-flavor (charm) jets [98]. Scale factors are applied to the simulated events to correct for differences in -tagging efficiency between data and simulation. These scale factors are measured as a function of jet T using a likelihood-based method in a sample highly enriched in¯events [98]. Only jets with | | < 2.5 are considered for flavor-tagging.
The energy of -tagged jets is corrected for the possible contribution of muons from semileptonic -hadron decays. In addition, the undetected energy of the neutrinos and out-of-cone effects are corrected for with scale factors derived as a function of the -jet T from a¯MC sample. The two corrections together improve the resolution of the invariant mass of the two jets with the highest -tagging score (¯) by about 20%. The procedure closely follows the one used in Ref. [102].
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 Δ = 0.4 of a selected photon. Next, jets within Δ = 0.2 of electrons are removed. In the last step, electrons and muons within Δ = 0.4 of any jet are removed.
The missing transverse momentum miss T in an event is calculated as the magnitude of the negative vectorial sum of the transverse momenta of all selected and calibrated physics objects that can be matched to the primary vertex. A component called the "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 miss T calculation [103].

Event selection
Both the nonresonant and resonant searches employ multivariate analysis techniques to select events. Events are selected if they satisfy a common set of preselection requirements; they are then required to fulfill different requirements for the nonresonant search and the resonant search. Both analyses employ a fit of the diphoton invariant mass distribution to extract the signal contribution.

Common preselection
For both the nonresonant and resonant analyses, events are selected using diphoton triggers requiring two reconstructed photon candidates with minimum transverse energies of 35 GeV for the leading photon and 25 GeV for the subleading photon, where leading (subleading) refers to the photon candidate with the highest (second-highest) transverse energy [104]. The triggers used in 2015 and 2016 required both photons to satisfy a Loose photon identification criterion, while the Medium criterion [90] was adopted in 2017-2018 to cope with the increased interaction rate. Once the full diphoton event selection described in this Section is applied, the average trigger efficiency for → events is found to be greater than 99% for the 2015-2016 data-taking period, and greater than 98% for the 2017-2018 data-taking period.
On top of the trigger requirements, events are selected if: • At least two photons satisfy the object selection criteria detailed in Section 4.1.
• The leading (subleading) photon T is larger than 35% (25%) of the mass of the diphoton system.
• Exactly two -tagged jets are present. In order to remain statistically independent of the ATLAS search for →¯¯ [105], any event with more than two -jets passing the 77% efficient working point is rejected.
• No electrons or muons are present.
• Fewer than six central (| | < 2.5) jets are present. This helps to reject¯events where the top quarks decay hadronically.
The acceptance times efficiency of the common preselection for the SM ggF simulation sample is 14% and it is 8.5% (14%) for a resonant scalar particle with = 300 GeV ( = 500 GeV). Multivariate techniques are then used to target the nonresonant ggF production mode or the resonant production mode of Higgs boson pairs. For both the nonresonant and resonant analyses, the resulting categories, each based on a boosted decision tree (BDT), must have at least nine continuum background events in the data sideband region. This guarantees that the final selection on the data sideband retains enough events to perform a fit to the distribution of the diphoton invariant mass, . The data sideband region is defined as the range 105 < < 160 GeV, excluding 120 < < 130 GeV.
The invariant mass of the diphoton plus b-tagged jets system, * ¯, is defined as * ¯=¯− − + 250 GeV (where 250 GeV is about twice the Higgs boson mass value and¯is the invariant mass of the two jets with the highest -tagging score). It is used to implement selection criteria for both the nonresonant and resonant analyses. Figure 4 shows that, compared with¯, the * ¯v ariable improves the four-object mass resolution, particularly for resonant signal particles decaying into a pair of Higgs bosons, due to detector resolution effects canceling out.

Nonresonant selection
Following the preselection, events are divided into two regions using the value of the * ¯v ariable. A high mass region, with * ¯> 350 GeV, targets the SM signal ( = 1), while a low mass region, with * ¯< 350 GeV, is used to retain sensitivity for BSM signals ( = 10). The dependence of * ¯o n can be seen in Figure 5. In each mass region, a dedicated BDT is trained using XGBoost [106] to discriminate between a benchmark signal and a combination of ,¯, , and simulated backgrounds. In the high mass region, the SM sample is used as signal, while in the low mass region, the = 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 -tagging information. Observables based on the kinematic properties of the reconstructed photons, such as the leading and subleading photon's angular information, and the transverse momentum of the diphoton system divided by its invariant mass, are combined with jet-based information. The "single topness" variable ( ) is also used. It is defined as = min where the minimum is taken over all combinations of three jets in the event (with no requirements on -tagging status), = 80 GeV, and = 173 GeV. Among the input variables in Table 2,¯and T show the highest discriminating power against the +jets continuum background. Particular care was taken to ensure that the BDT event selection does not lead to biases in the background distribution. Variables which have a strong correlation with the diphoton invariant mass are avoided in the training in order to prevent the BDT event selection from biasing the background distribution. To this end, the transverse momentum values of the photons are divided by before being used as BDT input variables. A check for potential biases in the background distribution is described in Section 5.2. The BDT score distributions in the low mass and high mass regions are shown in Figure 6 for events passing the common preselection. In each mass region, two categories based on the BDT score are defined. The boundaries of the categories are chosen by maximizing the combined number-counting significance [107] using signal and background yields in the mass window 120 < < 130 GeV in the chosen categories. The four resulting BDT categories are defined in Table 3. Table 2: Variables used in the BDT for the nonresonant analysis. All vectors in the event are rotated so that the leading photon is equal to zero, while their relative azimuthal angular differences are kept unchanged.

Resonant selection
The resonant search uses a multivariate analysis based on a BDT technique. A potential limitation of a BDT-based selection is the low number of background events for higher resonance masses. To overcome this limitation, BDTs are trained jointly for all resonance masses hypotheses and background. The signal events corresponding to the different mass hypotheses are combined. The signal events are reweighted event-by-event to match the * ¯d istribution of the background events, such that the training is independent of the resonant signal mass hypothesis, but it still retains information of the correlation with the rest of the event variables. The procedure allows to reduce the fluctuations of the BDT performance between nearby mass points. It was checked that the procedure provides similar or better performance than employing separated BDT trainings for each individual signal mass hypothesis.
Using the TMVA toolkit [108], two BDTs are trained to better separate the signal from backgrounds of different nature: the and¯backgrounds (BDT ) and the single-Higgs-boson background (BDT Single ), where the and¯processes produce the dominant resonant backgrounds. The nonresonant process is not included in the training of the BDTs. A complete list of the variables used for the BDT training is given in Table 4. The miss T information is used in the training because it is useful in rejecting the single-Higgs-boson background (¯in particular) and the¯background.

Variable Definition
Photon-related kinematic variables T , Transverse momentum and rapidity of the diphoton system

Missing transverse momentum
The combined BDT score of an event is obtained by combining the two BDT scores in quadrature, as shown in Eq. (2): The two coefficients 1 and 2 ( 2 = 1 − 1 ) and BDT tot take values in the range [0,1]. Only events passing a minimum requirement on the value of BDT tot are considered in the analysis. To reduce the effect of limited statistics of the Monte Carlo samples in the determination the minimum value required for BDT tot , and of the 1 and 2 coefficients, a two-stage procedure is used. The values of the 1 and 1 coefficient is sought (resulting in 1 = 0.65) across all the resonances so that the selection will have common coefficients for all resonance mass points, but different BDT tot values. For each of the resonance mass hypotheses, a requirement is set on the * ¯v alue to select events within ±2 of the expected mean value for signal events (where is the standard deviation parameter of the Crystal Ball function that best fits the * ¯d istribution). In the case of the 900 GeV and 1000 GeV mass hypotheses, the requirement is relaxed to ±4 to increase the number of events used for the signal extraction. The BDT tot distributions are shown in Figure 7 for two signal hypotheses ( = 300 GeV and 500 GeV).

Comparison of data and predictions
The analysis selection described in Section 4.2 requires two "tight" photons and this region is mainly composed of , -jet, and dĳet events, where either one or two jets are misidentified as photons. The fractional contribution of each component can be estimated using a data-driven approach [109] based on the photon identification and isolation distributions from genuine and misidentified photons. After the common preselection, (85 ± 3)% of sideband events are genuine diphoton events, with the remaining (15 ± 4)% consisting of -jet events and a negligible number of dĳet events. The uncertainties in the above fractions include both the statistical and systematic uncertainties, where the systematic uncertainty is estimated by using different photon identification criteria. In the BDT-based categories, the proportion of genuine diphoton events increases but the method suffers from a low event count for both the nonresonant and resonant cases. Figure 8 shows the level of agreement between data and the background prediction for the and * d istributions, after the common preselection. The continuum background is scaled by the , -jet, and dĳet fractions and normalized to the data sideband. The +jets continuum background is further divided according to the flavors of the two jets (for example¯or other jets). This decomposition is taken directly from the proportions predicted by the S event generator, as described in Section 3, and it is shown for illustration. Figures 9 and 10 show the distribution after the nonresonant and resonant BDT categorization and for two benchmark mass points = 300 GeV and = 500 GeV. The figures illustrate the signal and background composition. The background distributions shown in these figures are not directly used to model the background processes in the analysis workflow explained in Section 5.

Signal and background parameterization
The signal and backgrounds are extracted by fitting analytic functions to the diphoton invariant mass distribution in the range 105 < < 160 GeV in both the resonant and nonresonant searches.

→ parameterization
For the di-Higgs signal and single-Higgs-boson background distributions, 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 diphoton invariant mass distribution shapes are modeled with a double-sided Crystal Ball function [91, 110], which is characterized by a Gaussian core and asymmetric power-law tails. This function allows the modeling of event distributions in which non-Gaussian tails can arise from experimental effects, such as photon energy mismeasurements.
The shape parameters are determined by fitting the diphoton 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 nonresonant search, the parameterized form of is obtained from the simulation of the ggF and VBF processes with = 1, described in Section 3. No significant dependence of the functional form on was found. For the resonant search, the functional form is obtained from the simulation of the heavy-resonance signals. Table 5 shows   = 500 GeV for the resonant search. The data-derived fractions of nonresonant , -jet or jet-, and dĳet background are applied and the total background is normalized to the data sideband. The scalar resonance signal is scaled to a total production cross section ( → → ) = 370 fb for = 300 GeV or ( → → ) = 67 fb for = 500 GeV, corresponding to the expected exclusion limits presented in Section 7.3. the resolution parameter of the double-sided Crystal Ball functional form fit to the distribution for simulated Higgs boson pair events in the nonresonant 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 for both the single Higgs and di-Higgs boson processes.

Background parameterization
The continuum diphoton background is modeled using a functional form chosen by fitting a highly populated MC background template. Given the high purity quoted in Section 4.3, the background template is constructed in each category from the +jets simulation which is normalized to the data sideband in the mass windows of 105-120 GeV and 130-160 GeV in the spectrum.
The potential bias associated with the choice of a specific analytic function to model the continuum background is estimated for each category as prescribed in Refs. [91,111]. It is obtained as the signal event yield extracted from a signal-plus-background fit to the background-only diphoton invariant mass distribution in the range 105 < < 160 GeV. This bias is also called the "spurious signal". The number of fitted signal events is computed for Higgs boson masses in intervals of 1 GeV from 121 GeV to 129 GeV. The bias is taken as the largest number of fitted spurious 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 from the functions with a spurious signal smaller than 20% of the data's statistical uncertainty plus two times the MC statistical uncertainty. In each category of the nonresonant search and in all the analysis regions defined in the resonant search, an exponential function exp( · ) is found to be the best choice since it has the smallest number of degrees of freedom and yields a consistently small bias. Wald tests [112] on data show that the data do not prefer a higher-degree functional form to the exponential form. The bias is found to be at most 40% (60%) of the expected error on the fitted signal yield originating from the data statistics for the nonresonant (resonant) search. The difference in shape between the simulated events and the exponential form measured in the sidebands is found to have a negligible impact on the spurious signal. To study the impact on the spurious signal from the difference in shape between the MC background template and the data sideband, the MC template is reweighted to data using a linear function of derived from the MC-data difference. The spurious signal value evaluated from the reweighted MC template is found to be compatible with the nominal one used as the systematic uncertainty.

Systematic uncertainties
The sensitivity of the analysis is limited by the statistical precision. The assessment of the systematic uncertainties is described below and their impact on the results is discussed in Section 7.4. The uncertainty in the integrated luminosity of the full Run 2 data set is 1.7% [39], obtained using the LUCID-2 detector [38] for the primary luminosity measurements.
The continuum 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 propagated through the full analysis procedure, as described in the following.
The efficiency of the diphoton trigger used to select events is evaluated in simulation and data using a bootstrap method and radiative decays [113]. The difference between data and MC is taken as a systematic uncertainty. In the diphoton invariant mass range of 105 < < 160 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 → + − events in data with that in simulation [114]. The resulting uncertainty is found to have a negligible effect on the signal selection efficiency.
The systematic uncertainties from the photon identification and isolation efficiencies are estimated following the prescriptions in Ref. [90]. They are evaluated by varying the correction factors for photon selection efficiencies in simulation by the corresponding uncertainties and affect the diphoton selection efficiency. The experimental uncertainties in the photon energy scale and resolution are obtained from Ref. [90].
The jet energy scale and resolution uncertainties [115] affect the¯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 miss T calculation. In addition, the uncertainties in the scale and resolution of the miss T soft term are evaluated by using the method described in Ref. [103]. The flavor-tagging uncertainties for -and -jets are estimated in¯events [98,116], while the misidentification uncertainty of light-flavor jets is determined using dĳet events [117]. Additional uncertainties from the -jet momentum correction accounting for the presence of muons and neutrinos are found to be negligible.

For single Higgs boson and SM
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 SM ggF production arising from the choice of renormalization scheme and scale of the top-quark mass and their combination with those from factorization and renormalization scale variations are based on Ref. [17]. The uncertainties in the cross sections, including effects of PDF+ s uncertainties, and the uncertainties in the → and →¯branching fractions, are taken from Refs. [14,56]. The uncertainty in the value of the Higgs boson mass is considered [13]. An additional 100% uncertainty is assigned to the single-Higgs-boson ggF and VBF production modes and to Higgs boson production in association with a boson. This is motivated by studies of heavy-flavor production in association with top-quark pairs [118, 119] and of boson production in association with -jets [120]. No additional heavy-flavor uncertainty is assigned to the single-Higgs-boson¯and production modes, where the dominant heavy-flavor production is already accounted for in the LO process. For all samples, the uncertainty related to the choice of parton showering model is evaluated by comparing the predictions of the nominal MC samples and alternative samples using a different parton showering model. In addition, for the nonresonant production processes, a systematic uncertainty is assigned to the reweighting, as discussed in Section 3.
In the resonant search, the systematic uncertainty sources considered are the same as for the nonresonant case, except that uncertainties due to the finite order of the QCD calculations are not assigned to the resonant signal cross section. In this search the SM production processes are considered as background.

Results
The statistical framework used to derive the results for both the nonresonant and resonant searches is described in the following.

Statistical framework
For both the nonresonant and resonant searches, the results of the analysis are obtained from a maximumlikelihood fit of the distribution in the range 105 < < 160 GeV, performed simultaneously over all relevant categories described in Section 4.2. The likelihood function is defined in Eq. (3): where for each event in a category , is the observed number of events, is the expected number of events, is the value of the probability density function (pdf), are nuisance parameters, and ( ) are constraint pdfs for the nuisance parameters.
The expected number of events , defined in Eq. (4), is the sum of the expected yields from di-Higgs boson production processes, single Higgs boson production, the nonresonant background, as well as the spurious-signal uncertainty: where is the signal strength, SS,c represent the nuisance parameters associated with the background function bias and yield represent the nuisance parameters affecting the event 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 nonresonant search, res bkg,c = ,c , while in the case of the resonant analysis res bkg,c = ,c + SM ,c .
The probability density function represents the shape information. The sum of the double-sided Crystal Ball functions modeling the production processes, single Higgs boson production, and the spurious signal, and of the analytic function modeling the nonresonant background as described in Section 5.2, is shown in Eq. (5): where shape represent 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 initially set to values from simulation. The likelihood function includes all the nuisance parameters that describe the systematic uncertainties. The signal strength is a free parameter in the fit. The measurement of the parameter of interest is carried out using a statistical test based on the profile likelihood ratio [107], as shown in Eq. (6): where and are respectively the parameter of interest and the nuisance parameters. In the numerator, the nuisance parameters are set to their profiled valuesˆ( ), which maximize the likelihood function for a fixed value of the parameter of interest . In the denominator, the parameter of interest and the nuisance parameters are set to the valuesˆandˆ, respectively, which jointly maximize the likelihood.
In the absence of signal, exclusion limits are set on Higgs boson pair production in the¯final state. The limits for both nonresonant and resonant production are calculated using the CL S method [121], with the profile-likelihood-ratio-based test statistic˜, defined in Eq. (7) [107]: The asymptotic approximation [107] is used for the test-statistic distribution. The inaccuracy of the asymptotic approximation is checked with pseudo-experiment studies for the key results reported in this paper. Figure 11 shows the background fits to the data. No significant excess over the SM background expectations is found, as summarized in Table 6. The statistical analysis sets a 95% CL upper limit on the nonresonant 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.2 (5.7) times the SM prediction is obtained. The expected constraints are obtained for a background hypothesis excluding → production. 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. A check of the results that quantifies the upper limits by using pseudo-experiments is performed and the increase of the limit value relative to the asymptotic approximation is found to be less than 8% for both the observed and expected upper limits. A check of the expected upper limits using pre-fit values for the nuisance parameters associated with the systematic uncertainties is performed and an increase of 4% relative to the nominal result is found. The difference is dominated by the contribution from the spurious signal.

Nonresonant search results
Upper limits on the production cross section are also computed as a function of , as shown in Figure 12. For this purpose, single-Higgs-boson production cross sections and Higgs boson decay branching ratios are assumed to have SM values, and the coupling strength between the Higgs boson and other particles are also set to their SM values [55]. The theory uncertainties related to the signal cross section are not included.
The expected constraints on at 95% CL, as obtained with a background hypothesis excluding → production, are [−2.4, 7.7], whereas the observed constraints are [−1.5, 6.7] at 95% CL. The inclusion of the VBF production mode tightens the constraints by about 5% relative to an alternative fit considering only the ggF production mode. Table 6: The number of data events observed in the 120 < < 130 GeV window, the number of signal events expected for = 1 and for = 10, and events expected for single Higgs boson production (estimated using MC simulation), as well as for continuum background. For the single Higgs boson, "Rest" includes the following production modes: VBF, , , and . The values are obtained from a fit of the Asimov data set [107] generated under the SM signal-plus-background hypothesis, = 1. The continuum background component of the Asimov data set is obtained from the fit of the data sideband. The uncertainties in signals and single Higgs boson background include the systematic uncertainties discussed in Section 6. The uncertainty in the continuum background is given by the sum in quadrature of the statistical uncertainty from the fit to the data and the spurious-signal uncertainty. An alternative statistical analysis consists in determining the best-fit value of the coupling modifier. The best-fit value of and its uncertainty are obtained by means of a negative log-likelihood scan. The coupling strengths of the Higgs boson to fermions and gauge bosons are set to their SM values. The values of the negative log-likelihood ratio, −2 ln Λ( ), as a function of are shown in Figure 13. The Asimov data set [107] is generated under the SM signal-plus-background hypothesis, = 1. All systematic uncertainties, including those of the theoretical prediction of the production cross section, are included. The best-fit value corresponds to = 2.8 +2.0 −2.2 ( +3.8 −4.3 ) for the 1 (2 ) confidence interval. The expected value corresponds to = 1.0 +5.5 −2.4 ( +7.3 −4.2 ) for the 1 (2 ) confidence interval. The second minimum in the expected likelihood scan curve corresponds to a similar fitted signal yield with respect to the point at the first minimum, which is a consequence of a higher cross section, but lower acceptance and worse signal-to-background separation. The distribution has a different shape at each of the two minima, as shown in Figure 5. ATLAS √ s = 13 TeV, 139 fb −1 HH→bbγγ Observed limit (95% CL) Expected limit (95% CL) Expected limit ±1σ Expected limit ±2σ Theory prediction SM prediction   Figure 14 shows the fit to the data of the resonant search for two benchmark values of the mass of a hypothetical scalar particle. No significant excess over the SM background expectations is found, as shown in Table 7. Figure 15 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 640-44 fb (391-46 fb) in the range 251 ≤ ≤ 1000 GeV. A check on the upper limits using pseudo-expeirments is performed. For the expected limits, the results based on the pseudo-experiments are found to be up to 10% higher compared with those derived based on the asymptotic approximation. As for the observed limits, the pseudo-experiments yield typically 10% higher results compared with the asymptotic approximation in most of the range explored, and the difference increases to 15% for the = 700 GeV signal hypothesis. Table 7: The number of events observed in the 120 < < 130 GeV window in data, the number of events expected for scalar resonance signals of masses = 300 GeV and = 500 GeV assuming a total production cross section ( → → ) equal to the observed exclusion limits of Figure 15, and events expected for SM and single Higgs boson production (estimated using MC simulation), as well as for continuum background. The values are obtained from a fit of the Asimov data set [107] generated under the signal-plus-background hypothesis. The continuum background component of the Asimov data set is obtained from the fit of the data sideband. The uncertainties in the resonant signals and the SM and single-Higgs-boson backgrounds include the systematic uncertainties discussed in Section 6. The uncertainty in the continuum background is given by the sum in quadrature of the statistical uncertainty from the fit to the data and the spurious-signal uncertainty. Expected limit (95% CL) Observed limit (95% CL) σ 1 ± Expected limit σ 2 ± Expected limit Figure 15: Observed and expected limits at 95% CL on the production cross section of a narrow-width scalar resonance as a function of the mass of the hypothetical scalar particle. The black solid line represents the observed upper limits. The dashed line represents the expected upper limits. The ±1 and ±2 variations about the expected limit due to statistical and systematic uncertainties are also shown.

Impact of systematic uncertainties
The dominant systematic uncertainties are listed in Table 8 for both the nonresonant and resonant searches. The main uncertainties are related to the choice of functional form for the continuum background (spurious signal), to the parton showering model, and to the photon energy resolution. Table 8: Breakdown of the dominant systematic uncertainties. The impact of the uncertainties corresponds to the relative variation of the expected upper limit on the cross section when re-evaluating the profile likelihood ratio after fixing the nuisance parameter in question to its best-fit value, while all remaining nuisance parameters remain free to float. The impact is shown in %. Only systematic uncertainties with an impact of at least 0.2% are shown. Uncertainties of the "Norm. + Shape" type affect both the normalization and the parameters of the functional form. The rest of the uncertainties affect only the yields.

Conclusions
Searches for nonresonant and resonant Higgs boson pair production are performed in the¯final state using 139 fb −1 of 13 TeV collision data collected with the ATLAS detector at the LHC. No significant excess above the Standard Model background expectation is observed. A 95% CL upper limit of 130 fb is set on the → nonresonant production cross section, where the expected limit is 180 fb. The observed (expected) limit corresponds to 4.2 (5.7) times the cross section predicted by the Standard Model. Constraints on the Higgs boson self-coupling are also derived and limits of −1.5 < < 6.7 are obtained, where −2.4 < < 7.7 is expected. The expected constraints on the nonresonant production cross section and on are obtained with a background hypothesis excluding → production. For resonant production of a scalar particle → →¯, upper limits on the production cross section are obtained for the narrow-width hypothesis as a function of . The observed (expected) upper limits are in the range 640-44 fb (391-46 fb) for 251 ≤ ≤ 1000 GeV. Compared to the previous ATLAS result based on 36 fb −1 of 13 TeV collisions, the present analysis uses a data set about four times larger, incorporates a categorization based on * ¯a nd multivariate event selections, and expands the analyzed mass range of the resonance search to lower values. The results improve upon the previous ATLAS limits on the →¯production cross section by up to a factor of five, and the allowed range shrinks by about a factor of two. For the resonant search, the expected limit on the cross section improves by a factor of two to three depending on the value. Of those improvements, a factor of two arises from the increase in integrated luminosity, while the additional improvement can be attributed to the use of multivariate techniques, more precise object reconstruction and calibration and, for the nonresonant search, the categorization based on * ¯. The present analysis also sets constraints that are tighter than those from a combination of ATLAS searches for production in up to 36 fb −1 of 13 TeV data.
[27] L. Randall     [93] ATLAS Collaboration, Muon reconstruction and identification efficiency in ATLAS using the full Run 2 collision data set at