Studies of new Higgs boson interactions through nonresonant $HH$ production in the $b\bar{b}\gamma\gamma$ final state in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector

A search for nonresonant Higgs boson pair production in the $b\bar{b}\gamma\gamma$ final state is performed using 140 fb$^{-1}$ of proton-proton collisions at a centre-of-mass energy of 13 TeV recorded by the ATLAS detector at the CERN Large Hadron Collider. This analysis supersedes and expands upon the previous nonresonant ATLAS results in this final state based on the same data sample. The analysis strategy is optimised to probe anomalous values not only of the Higgs ($H$) boson self-coupling modifier $\kappa_\lambda$ but also of the quartic $HHVV$ ($V=W,Z$) coupling modifier $\kappa_{2V}$. No significant excess above the expected background from Standard Model processes is observed. An observed upper limit $\mu_{HH}<4.0$ is set at 95% confidence level on the Higgs boson pair production cross-section normalised to its Standard Model prediction. The 95% confidence intervals for the coupling modifiers are $-1.4<\kappa_\lambda<6.9$ and $-0.5<\kappa_{2V}<2.7$, assuming all other Higgs boson couplings except the one under study are fixed to the Standard Model predictions. The results are interpreted in the Standard Model effective field theory and Higgs effective field theory frameworks in terms of constraints on the couplings of anomalous Higgs boson (self-)interactions.

The current measurements provide constraints on the strengths of the couplings of the Higgs boson to the heaviest of the SM elementary particles, and on the Higgs boson mass   .The latter is one of the parameters of the Higgs boson potential  () = 1  2  2   2 +      3 + 1 4       4 , where  ≈ 246 GeV is the vacuum expectation value of the Higgs field.Among the SM predictions of the Higgs sector that still remain to be verified are those for the coupling strengths of the interactions involving multiple Higgs bosons, such as the trilinear and quartic Higgs boson self-couplings,     and      , as well as the quartic couplings between two Higgs bosons and two  or  bosons,    ( = , ).In the SM, the trilinear and quartic self-couplings have the value  SM    =  SM     =  2  /2 2 , while the couplings    are related to the  and   couplings   through  SM   =  SM  /2.A significant effort has been dedicated by the ATLAS and CMS collaborations to search for processes that are particularly sensitive to     and    , such as Higgs boson pair production in gluon-gluon fusion (ggF) and vector-boson fusion (VBF).In the SM, ggF  production proceeds through the destructive interference of two leading Feynman diagrams: one for the process  →  * → , involving an intermediate virtual Higgs boson ( * ) and a  vertex (Figure 1(a)), and a second one describing a loop-mediated process in which two Higgs bosons are radiated off a virtual quark (Figure 1(b)).VBF  production is induced at tree level in the SM by three Feynman diagrams in which the two vector bosons radiated by the scattering quarks either fuse into a virtual Higgs boson  * decaying into two Higgs bosons via a  vertex (Figure 2(a)), fuse into two Higgs bosons via a  vertex (Figure 2(b)), or produce two Higgs bosons via -channel scattering through two  interactions (Figure 2(c)).The amplitudes of diagrams involving a  vertex are proportional to     , while those of diagrams involving a  vertex are proportional to    .For this reason, the results of the searches for  production can be used to infer the values of the coupling modifiers   ≡     / SM    and  2 =    / SM   .An observed value of these coupling modifiers significantly different from unity would provide a proof of non-SM Higgs boson interactions [20].q q H V V κ λ κ V (a) Trilinear coupling In the SM, these processes are expected to be rare, with cross-sections that are about three orders of magnitude smaller than those of single Higgs boson production:    ggF = 31.1 +2.1 −7.2 fb [21][22][23][24][25][26][27][28] and    VBF = 1.73 ± 0.04 fb [29][30][31] for   = 125 GeV and a proton-proton centre-of-mass energy of √  = 13 TeV.It is thus crucial to analyse the latest available data sample and reconstruct as many decay final states of the Higgs boson pairs as possible.The most stringent constraints on   and  2 exploit the entire sample of proton-proton ( ) collisions provided by the LHC during its second phase of data-taking (Run 2, 2015-2018), and a multitude of Higgs boson decay channels.In particular, the ATLAS experiment recently released the results of searches based on the full Run 2 data in the three most sensitive channels,  b [32],  b +  − [33], and  b b [34], and their combination [35].No excess over the SM background was observed, and constraints on the coupling modifiers   and  2 were set at the 95% confidence level (CL).The observed (expected for   = 1) 95% confidence interval for   when all other coupling strength modifiers are set to unity is −0.6 <   < 6.6 (−2.1 <   < 7.8) after combining the three  decay channels.For  2 , the observed (expected) 95% confidence interval when all other coupling strength modifiers are set to unity is 0.1 <  2 < 2.0 (0.0 <  2 < 2.1).With a similar data sample, CMS also reported similar results in their  b [36],  b +  − [37], and  b b [38,39] channels, observing 95% CL intervals of −1.2 <   < 6.5 and 0.7 <  2 < 1.4, based on a different statistical procedure [8].
The  b final state has an expected branching ratio (0.26%) that is significantly smaller than that of  b b (34%) and  b +  − (7.3%).However, the larger expected signal-to-background (S/B) ratio and the higher trigger efficiency and thus larger acceptance in phase-space regions (e.g., at small values of the  invariant mass), where potential deviations from the SM are expected to be enhanced, compensate for the lower expected event yield and lead to a sensitivity similar to that of the other two decay modes.The latest results for  →  b with the full Run 2 ATLAS data, based on the event selection and classification of Ref. [32] but using the statistical procedures of Ref. [35] and of this analysis, yields the following observed (expected) one-dimensional 95% confidence intervals: −1.4 <   < 6.5 (−3.2 <   < 8.1) and −0.8 <  2 < 3.0 (−1.6 <  2 < 3.7).This paper presents an updated search for nonresonant Higgs boson pair production in the  b final state using the full Run 2 ATLAS data, superseding and expanding upon the nonresonant results of Ref. [32].Compared to the previous publication, an identical event selection and a similar analysis strategy are used, but a reoptimised classification of events in categories with different S/B leads to a higher sensitivity to the   and  2 coupling modifiers.The new event classification relies on improved multivariate classifiers, also exploiting the kinematic features of VBF  production for SM and anomalous values of   and  2 .After the events are classified in mutually orthogonal event categories, the signal cross-section is estimated through a simultaneous maximum-likelihood fit to the diphoton invariant mass spectrum of the selected events in each category.The fit probes an enhancement in event yields around the experimental value of the Higgs boson mass over the predicted background, consisting of the sum of a monotonically decreasing distribution from continuum photon and jet production and a peak from singly produced Higgs bosons decaying into two photons.
Another novelty compared to the previous publication is the interpretation of the results in two effective field theory (EFT) extensions to the SM, the Higgs effective field theory (HEFT) [40,41] and the SM effective field theory (SMEFT) [42,43].The data are used to set constraints on the Wilson coefficients of operators of the EFT Lagrangians describing anomalous Higgs boson interactions in both frameworks.HEFT and SMEFT describe the same effective interactions, but with different bases of operators.One advantage of HEFT compared with SMEFT is that it provides a one-to-one relation between operators (and corresponding Wilson coefficients) and effective interactions, which allows single-and di-Higgs boson couplings to be separated, leading to simplified  interpretations.This paper is organised as follows.Section 2 describes the experimental apparatus.The data and simulated event samples used for the measurements are summarised in Section 3. Section 4 is devoted to the event selection and classification, with an emphasis on the novelties of the latter compared to the previous publication.The   signal and background models used in the final fit are described in Section 5.The systematic uncertainties in the measurement and the results are given in Sections 6 and 7, respectively.Finally, the procedure and results of the EFT interpretation are detailed in Section 8. Section 9 provides the conclusions.

The ATLAS detector
The ATLAS detector [44] 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 of || < 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 [45,46].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 || = 2.0.The TRT also provides electron identification information based on the fraction of hits (typically 30 in total) above a higher energy-deposit threshold corresponding to transition radiation.
The calorimeter system covers the pseudorapidity range of || < 4.9.In 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 with || < 1.7, and two copper/LAr hadron endcap calorimeters.The solid angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules optimised 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.Three layers of precision chambers, each consisting of layers of monitored drift tubes, cover the region || < 2.7, complemented by cathode-strip chambers in the forward region, where the background is highest.The muon trigger system covers the range of || < 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 [47].The first-level 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 to disk at about 1 kHz.
An extensive software suite [48] is used in data simulation, in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment. 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe.The -axis points from the IP to the centre 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 :  =

Data and simulation samples
The measurements presented in this paper use   collision data collected by the ATLAS experiment during the LHC Run 2 at √  = 13 TeV.After data quality requirements [49], the integrated luminosity of the data sample is 140.1 ± 1.2 fb −1 [50].
The simulated event samples used in this study are summarised in Table 1.Besides the samples already used in Ref. [32], VBF  samples were produced for additional  2 and   variations (where   =   / SM  is the  coupling modifier), and a dedicated diphoton + two -jet sample was generated.
Table 1: Summary of the nominal Higgs boson pair signal, single Higgs boson background and continuum background event samples used in this analysis.The generator used in the simulation, the parton distribution function (PDF) set, and the set of tuned parameters (tune) are also provided.The final two columns list the accuracy in QCD of the event generator and the order in QCD of the calculated cross-section for the  signal and the single Higgs boson background (LO: leading order, NLO: next-to-leading order, NNLO: next-to-next-to-leading order, N 3 LO: next-to-next-to-next-to-leading order).More details are given in the text and in Ref. [32].The accuracy and cross-sections for the nonresonant background processes are omitted since their shape parameters and overall normalisation are determined from fits to the data.[69].The particle-level    spectrum for any generic value of   is calculated from the    distributions of three ggF  samples generated at particle level for   = 0, 1, and 20.To determine the potential 'non-closure' in the reweighting process from residual kinematic effects, the procedure is validated by comparing the predicted event yields and kinematic distributions of the simulated sample generated with   = 1 and reweighted to   = 10 with those of the simulated sample generated under the hypothesis   = 10.Furthermore, 12 additional VBF  samples were generated and simulated with the same set-up and settings as the nominal VBF sample but using non-SM combinations of the coupling strength scale factors   ,  2 and   .A linear combination of a 'basis' formed by the SM sample and five of the other 12 samples, corresponding to the combinations of the   ,  2 , and   couplings (1, 1.5, 1), (0, 1, 1), (10, 1, 1), (1, 3, 1), (−5, 1, 0.5), is used to determine the expected yields and distributions for any value of   ,  2 , and   .The remaining seven samples are compared with the corresponding predictions from the interpolation procedure for validation purposes.The same procedure was used in the measurements presented in Refs.[34,35].
Background samples include simulated events of single Higgs bosons decaying into  produced by ggF, VBF, in association with a  or  boson, with a  t or  b pair, or with a single top-quark .Simulated event samples of continuum diphoton production in association with top quark pairs ( t) or with jets from quarks of other flavours (+jets) were also produced, to optimise the event classification described in Section 4.2.In addition to the previous samples, shared with Ref. [32], a sample of simulated continuum diphoton plus two -jets events ( b) was generated with Sherpa 2.2.12 [68] using NLO matrix elements for the production of the two photons and the two -quarks in the four-flavour scheme, with additional jets produced in the parton shower.Due to the increased efficiency from generator-level requirements on the -quarks, the use of this new sample reduces the statistical uncertainty in the main component of the nonresonant background originating from  b events by a factor of two, despite containing about 60 times fewer simulated events than the inclusive diphoton sample.This sample was used to study the background diphoton invariant mass distribution, as described in Section 5.
All generated samples were passed through a detailed simulation of the ATLAS detector response [70] based on Geant4 [71], except for the inclusive diphoton sample, which was interfaced to a fast detector simulation based on a parametric description of the calorimeter response [72], and for the ggF  particle-level samples used for the    -based reweighting procedure, for which the detector response was not simulated.The generation of the simulated event samples includes the effect of multiple inelastic   interactions per bunch crossing, and the effect on the detector response of interactions from bunch crossings before or after the one containing the hard interaction.The inelastic   events were generated with Pythia 8.186 using the NNPDF2.3loPDF set and the A3 tune [73].The Higgs boson mass was assumed to be 125 GeV in both simulation and the analysis of the data.The impacts of the differences relative to the best-fit values of the   measurements reported in Refs.[18,19], and the effects of the corresponding experimental uncertainties in   , are negligible.

Event selection and classification
The same preselection as described in Ref. [32] is used to suppress the background while providing good signal efficiency.It is briefly summarised in Section 4.1.The selected events are then classified into orthogonal categories based on multivariate discriminants using several input kinematic quantities.The definition of the event categories, described in Section 4.2, is chosen in order to optimise the expected constraints on the coupling modifiers   and  2 .

Event preselection
To identify  →  decays, events were collected with diphoton triggers [74] with nominal transverse momentum ( T ) thresholds of 35 GeV and 25 GeV for the leading-and subleading- T candidates, respectively.Selected events are required to contain two photon candidates in the acceptance of the finely segmented part of the electromagnetic calorimeter (|| < 1.37 or 1.52 < || < 2.37).The candidates must be identified as photons by an algorithm based on the shower shapes reconstructed in the calorimeter.Of all potential reconstructed collision vertices, the primary diphoton vertex (PV) is selected by a neural-network algorithm using extrapolated photon trajectories and tracks associated with the candidate vertices [75].
The photon candidates must also meet the requirements of an isolation algorithm based on the energy flow in the calorimeter and the total transverse momentum of charged particle tracks from the PV in the inner detector, in cones surrounding the photon direction [76].The two leading photons passing these selections are then required to have an invariant mass   between 105 and 160 GeV and transverse momenta above 35% and 25% of   .
Jets are reconstructed from particle-flow objects built from noise-suppressed positive-energy topological clusters in the calorimeter and reconstructed tracks using the anti-  clustering algorithm with the parameter  = 0.4 [77,78].Jet candidates are required to have  T > 25 GeV and || < 4.4.Jets in the fiducial acceptance of the inner detector (|| < 2.4) and with  T < 60 GeV must be identified by a 'jet-vertex tagger' as originating from the PV [79].To target  →  b decays, events are required to contain exactly two -tagged jets, defined as central jets (those in the acceptance of the inner detector (|| < 2.5)) that satisfy the criteria of the 'DL1r' -tagging algorithm with a nominal efficiency of 77% for -jets and a misidentification rate of 1/170 (1/5) for light-flavour (charm) jets in  t simulated events [80].A correction factor is applied to the energy of the two -tagged jets to account for possible contributions from muons originating from semileptonic -hadron decays and undetected energy from neutrinos and out-of-cone effects [32].Jets failing to satisfy the -tagging requirement are ranked from first to last based on a discrete -tagging score defined by three bins, corresponding to central jets with DL1r efficiencies of 77%-85% and 85%-100%, and non-central jets.Jets with the same score are ranked by  T .
Events with six or more central jets, or with one or more isolated lepton (electron or muon) candidates with  T > 10 GeV and passing the lepton identification criteria are rejected in order to suppress background from  t () and inclusive  t production.No requirements are made on the number of non-central jets.
The efficiency of the event preselection is 13% (9%) for SM ggF (VBF)  events.The number of events selected in data in this inclusive signal region is 1874.With this selection, approximately 45% of the continuum background consists of events with two genuine -jets and two prompt photons, 40% consists of events with two genuine prompt photons and at least one misidentified -jet, and 15% consists of events with at least one misidentified photon.

Event categories
The kinematic properties of improves the signal mass resolution due to the cancellation of detector resolution effects [32].
In each of the two  *  b regions, a dedicated boosted-decision-tree (BDT) discriminant is trained to distinguish  signals from the background arising from  →  decays in single Higgs boson production events and from the continuum diphoton background from  t and +jets events.The training is performed with the XGBoost program [81] using only simulated event samples.In the high mass region, the signal samples used for training include SM ggF and VBF  events, as well as the five non-SM samples of the VBF  basis.In the low mass region, the signal samples consist of non-SM ggF  events corresponding to   = 10 and   = 5.6, plus the same five non-SM VBF  basis samples.The choice of   = 5.6 corresponds to a large anomalous value of   that is not yet excluded with a high confidence level by the previous search in this channel.However, it is observed that the training is relatively stable for variations of the order of unity on the   value used in training.
The BDT discriminant uses the same input variables that were used for the analogous multivariate discriminant in Ref. [32] (denoted by baseline variables), complemented by a set of additional observables that provide further discrimination between the background and the signal, mainly from VBF  production.The baseline variables include kinematic properties of the two photon and the two -jet candidates, the scalar sum  T of the  T of all the jets, and the magnitude  miss T and direction  miss of the missing transverse momentum vector ì  miss T [82].Another baseline variable is the single-topness   , quantifying how likely any three-jet combination in the event is to originate from a  →   →  q′  decay: where   and   are the masses of the  boson and of the top quark, and the minimum is evaluated over all combinations of any three jets in the event, with no requirements on whether they are -tagged.
The additional variables include, for events with at least four jets, the  T , , , and discrete -tagging score of the third and fourth jets.Events with at least four jets can arise from VBF  production, in which the scattered quarks responsible for the VBF process hadronise after having radiated a weak boson and produce two forward, high-momentum jets ('VBF jets').In events with exactly four selected jets, the two non -tagged jets are considered as VBF-jet candidates.In events with at least five selected jets (about 25% of the VBF  events passing the previous requirements according to the simulation), the two non -tagged jets that are considered as VBF-jet candidates are determined by means of a BDT classifier ('VBF-jet tagger').The inputs of the VBF-jet tagger consist of: (i) for each non -tagged jet , its  T , , and Δ and Δ separations from the  b system; (ii) for each   pair, its invariant mass, Δ between the two jets, Δ and Δ separations from the  b system, and  T , , and invariant mass of the  b   system.The BDT is trained on simulated SM VBF  events using the pair of jets matched to the scattered quarks as signal, and all other pairs of jets as background.After training, the VBF-jet tagger is applied to all possible jet pair combinations in data and simulated events, and the jets belonging to the pair with the highest tagger score are considered as VBF-jet candidates.Their invariant mass and pseudorapidity difference are then used as input variables for the event classification BDTs.In simulated VBF  events with at least three non -tagged jets, the VBF-jet tagger is able to correctly identify the VBF-jet pair in 95% of events.
A second set of additional variables used as input to the event classification BDTs consists of event-level kinematic quantities such as  *  b and the angular separation Δ(, ) (Δ(, b)) between the two photon (-tagged jet) candidates.Finally, three event-shape observables are also used: the transverse sphericity  ⊥ [83], the planar flow   [84], and the transverse momentum balance, defined as The relative weights of the training samples, as well as the values of the XGBoost hyperparameters, are tuned using a Bayesian optimisation algorithm that maximises the expected combined number-counting significance  [85] of a benchmark signal using the signal and background yields in each category in the diphoton invariant mass range 120 GeV <   < 130 GeV, as described below.
After training, three categories (labelled 'High Mass ',  = 1..3) in the high mass region and four categories (labelled 'Low Mass ',  = 1..4) in the low mass region are defined based on the high mass region and low mass region BDT discriminants, with a higher category index  corresponding to higher BDT scores and more signal-like events.Events from the inclusive signal region are thus classified in seven orthogonal exclusive signal regions based on the value of  *  b and of the BDT scores.Events with a BDT score lower than the threshold defining the category with the lowest index in the corresponding low or high mass region are discarded.The values of the BDT scores used to define the categories are chosen by maximising the combined number-counting significance of all categories in a region for a benchmark signal using expected signal and background yields in the diphoton invariant mass range 120 GeV <   < 130 GeV.During this optimisation process, each category must contain at least nine expected continuum background events in the   sidebands, i.e. excluding the region 120 GeV <   < 130 GeV, in order to have sufficient events to constrain the shape of the diphoton invariant mass distribution of the continuum background when the selection is evaluated on the data.In the high mass region, the signal yield is computed from the sum of the expected SM ggF and VBF  contributions, while in low mass region, the signal yield is computed from the ggF    = 5.6 and VBF    = 10 predictions.
The BDT discriminant distributions in the low and high mass regions observed in data in the   sidebands are shown in Figure 3. Also illustrated for comparison are the expected BDT score distributions for the dominant nonresonant background from the  + jets sample, the resonant single Higgs boson background, and the ggF and VBF  signals for different values of   and  2 .The values of the BDT scores that define the categories are represented by vertical dashed lines.In total, 340 events in the range of 105 GeV <   < 160 GeV are retained from the 1874 passing the initial preselection.

Signal and background modelling of the diphoton mass spectrum
The signal, resonant and nonresonant background yields in each category are determined from unbinned fits to the diphoton invariant mass distributions in the signal regions, as described in Section 7. The signal and background   distributions in each category are independently modelled by means of analytical functions chosen as follows.
The   distributions of signal events and resonant background events from single Higgs bosons decaying into  are described by double-sided Crystal Ball functions [75,86].The shape parameters are obtained from fits to simulated SM  events, and then either fixed in the final fits (parameters describing the tail of the distribution) or constrained around the initial values within the uncertainties resulting from the photon energy calibration.The same model is found to describe selected single Higgs boson and Higgs boson pair events well for both SM and non-SM coupling values.Signal + background fits performed on a combination of signal and resonant background events from simulation and the expected nonresonant background distribution show negligible signal yield non-closure resulting from this assumption.
The   distributions of the nonresonant diphoton background are modelled with exponential functions, whose normalisation and shape parameters are obtained from the fit to the data.The chosen exponential model in each category has two degrees of freedom and is found to describe the data well in the   Figure 3: BDT score distributions for simulated ggF and VBF  →  b signal events and simulated background events from nonresonant  + jets and singly produced Higgs bosons decaying into  for the (a) low and (b) high mass regions.The data in the   sidebands, which are not expected to be populated by single nor double Higgs boson events, are also shown compared with the  + jets sample.The latter comprises the majority of the nonresonant diphoton background and is used in the training of the BDT.All distributions are normalised to unity.The vertical dashed lines correspond to the thresholds used to define the event categories.Events with a BDT score between 0 and the lowest threshold (thick dashed line) are discarded.Events satisfying the lowest threshold are categorised as Low Mass ,  = 1..4 (High Mass ,  = 1..3), with a higher category index  corresponding to higher BDT scores and more signal-like events.sidebands, as well as the background-only template obtained with the Sherpa 2.2.12  b sample normalised to the data in the   sidebands.The spurious signal [75,87] is defined as the maximum absolute value of the bias on the fitted signal yield in multiple signal+background fits to the background-only template, performed with varying mass assumptions on   ∈ [123, 127] GeV in intervals of 0.5 GeV.For each of the exponential models, the spurious signal is smaller than 20% of the statistical uncertainty in the expected fitted signal yield, plus twice the statistical uncertainty in the spurious signal itself.Alternative models with the same numbers of degrees of freedom, such as power functions of   , performed similarly to the exponential model.While the nominal templates are constructed with simulated  b events only, alternative templates accounting for potential shape differences due to other background components such as  q ( ≠ ),  , and   do not significantly alter the spurious signal value or the quality of the exponential fit.

Systematic uncertainties
Systematic uncertainties affect the shape and normalisation of the diphoton invariant mass distributions of the Higgs boson pair signal and single Higgs boson backgrounds.Nevertheless, due to the limited number of events and the small signal-to-background ratio, the impact of the systematic uncertainties is small compared with that of the statistical uncertainties.
The systematic uncertainties are computed separately for the ggF and VBF  production modes and for the various single Higgs boson production modes.Those from the same source are correlated between processes.For the ggF (VBF)  signal, for each source of uncertainty the corresponding estimate is obtained by taking the envelope of values computed using both the SM and the   = 10 simulated event sample (using the six VBF basis simulated event samples).
The uncertainty in the full Run 2 integrated luminosity is derived from dedicated measurements [50] using the LUCID-2 [88] detector.The diphoton trigger efficiency uncertainty is evaluated using radiative  boson decays and with events collected using prescaled lower-threshold triggers [74].The uncertainty in the vertex selection efficiency is evaluated by comparing the reconstruction efficiency of photon-pointing vertices in  →  +  − events in data with that in simulation [89].
The uncertainties in photon identification and isolation efficiencies are determined from control samples of prompt photons from photon+jet production and from radiative  boson decays and electrons [76].The uncertainties in the photon energy scale and resolution are determined from control samples of electrons from  boson and / decays and of photons from radiative  boson decays [76].
The uncertainties in the jet energy scale and resolution are determined from control samples of jets recoiling against well calibrated particles such as photons,  bosons or already calibrated jets [90].Additional uncertainties from the simulation account for potential differences between the response for -jets and jets from gluons and light quarks.The uncertainties in the flavour-tagging efficiencies and misidentification rates are estimated by using  t events for -and -jets and +jets events for light-flavour jets [80,91,92].Theoretical uncertainties due to missing higher-order terms in the perturbative expansion of the crosssection, the PDF set, and the value of  s affect the total expected yields of single Higgs boson and Higgs boson pair events, and their fractional contributions to each category.These uncertainties are evaluated by considering alternative choices of factorisation and renormalisation scales, PDF sets, and the value of  s .For SM Higgs boson pair production, the values of the QCD scale and PDF+ s total cross-section uncertainties are taken from Ref. [93].For SM  production through ggF, the QCD scale and PDF+ s cross-section uncertainties are further combined with the top-quark mass scale uncertainty according to the prescription described in Ref. [28].The uncertainties in the  →  and  →  b branching ratios are also included [94].
For the signal and the ggF, VBF, and  t single Higgs boson processes, the uncertainty due to the choice of parton shower model is evaluated by comparing the predictions of the nominal simulation using the Pythia 8 model with an alternative simulation in which the same generator-level events are showered with Herwig 7.An additional 100% uncertainty in the yields of single Higgs boson ggF, VBF and   production modes is applied, motivated by studies of heavy-flavour production in association with top-quark pairs [95,96] and  boson production in association with -jets [97].
For the ggF  process with   ≠ 1, a systematic uncertainty is assigned to the   reweighting procedure by computing for each category the maximum deviation between the expected yields determined from the ggF  sample generated with   = 10 and the sample generated with   = 1 and reweighted to   = 10.For the VBF  process, a similar uncertainty for the potential non-closure of the procedure used to calculate the expected yield for any value of   and  2 from a linear combination of the six basis samples is determined for each category.It is calculated as the maximum difference, for the seven validation samples described in Section 3, between the expected yield calculated with the validation sample and that obtained from the linear combination approach.
An additional uncertainty in the signal yield is due to the choice of the background model and is assumed to be equal to the spurious signal described in Section 5.The larger equivalent integrated luminosity of the Sherpa 2.2.12  b sample used to create the background-only template, compared with that of the Sherpa 2.2.4 +jets sample used in the previous search published in Ref. [32], is more effective at suppressing statistical fluctuations in the template that would otherwise lead to overestimated spurious signals.As a consequence, the spurious signal obtained with the background template from the Sherpa 2.2.12  b sample in each category ranges between 10% and 50% of that from the background template produced with the Sherpa 2.2.4 +jets sample.Its impact on the expected upper limit on the  signal strength    , defined as the ratio of the Higgs boson pair production cross-section to its SM prediction, is thus at the permille level, compared to 3% in the previous analysis.
The impacts of the systematic uncertainties in the expected 95% CL upper limit on    , determined with the statistical interpretation described in the next section, are listed in Table 2.
Table 2: Breakdown of the dominant systematic uncertainties in the expected    upper limit at 95% CL.The impact of the uncertainties corresponds to the relative variation of the expected upper limit 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.Only systematic uncertainties with an impact of at least 0.1% are shown.

Systematic uncertainty source
Relative

Results
The results are derived using the statistical procedures outlined in Refs.[32,35,98] from the global likelihood function (, ).The set  contains the parameters of interest (POI) of the measurement, while  is the ensemble of nuisance parameters, corresponding to systematic uncertainties constrained by auxiliary measurements in control regions or by theoretical predictions, or to parameters such as the continuum background yields that are a priori unconstrained.The function (, ) is the product of the likelihood functions in each of the seven orthogonal categories, and of constraint terms for the nuisance parameters that are not freely floating in the fit.For each category, the likelihood function is determined from the corresponding signal and background models of the   probability density functions described in Section 5, the signal and background yield expectations for given values of  and , and the observed   distribution in data.
The constraints on the coupling strength parameters, expressed as 68% and 95% CL intervals, are determined with the same procedure as that of Ref. [35], using a profile-likelihood-ratio test statistic Λ(, ) computed from the likelihood function in the asymptotic approximation [85], where the POIs in  are the coupling strength modifiers .Signal strength upper limits are derived as in Ref. [32] using the CL s approach [99] from a separate test statistic q  that evaluates to zero when the parameter of interest  =    corresponding to the cross-section under study is lower than its maximum likelihood estimate (MLE) μ .The allowed   interval published in Ref. [32] was determined in a different way, from the range of   values for which the predicted  cross-sections are lower than the observed upper limits.The expected results are obtained with Asimov datasets [85] generated from the likelihood function after setting all nuisance parameters to their MLE in the fit to the data and fixing the POIs to the values corresponding to the hypothesis under test.The asymptotic results are found to agree within 10% with values obtained using pseudo-experiments.
Figure 4 shows the result of a background-only fit to the data, using the likelihood function  after fixing the parameters of interest corresponding to setting the signal cross-sections to zero.Table 3 compares the number of events in the observed data to the expected values in each category.No significant excess over the expected background is found, and a 95% CL upper limit of 4.0 on the total  production signal strength    (where only ggF and VBF processes are considered) is set, to be compared with an expected limit of 5.0 (6.4) in the background-only    = 0 (SM    = 1) hypothesis.If the VBF (ggF)  signal strength is fixed to the SM prediction, the observed upper limit on the ggF (VBF)  signal strength is 4.1 (96), while the expected upper limit, computed assuming  ggF = 0 ( VBF = 0), is 5.3 (145).The observed limits are tighter than the expected ones due to deficits in the signal regions of the most sensitive categories, as shown in Table 3.The compatibility between the best-fit value of    and the SM expectation is approximately 1.3 standard deviations.
The values of −2 ln Λ as a function of the coupling strength factor   or  2 under the hypothesis that all other coupling modifiers are equal to their SM predictions are shown in Figure 5.The observed (expected) constraints under this hypothesis are −1.4 <   < 6.9 (−2.8 <   < 7.8) and −0.5 <  2 < 2.7 (−1.1 <  2 < 3.3) at 95% CL.Two-dimensional constraints at 68% and 95% CL in the (  ,  2 ) plane are also shown in Figure 6, when all the other coupling modifiers are fixed to their SM predictions.
The impact of the systematic uncertainties on the results is small, leading to an increase of the upper limits on the signal strengths by 6%-7% and to a widening of 95% CL confidence intervals for the coupling modifiers by 2%-3% relative to the case in which systematic uncertainties are neglected.
Compared to the previous analysis of Ref. [32,35], the new event classification procedure leads to a reduction in the expected upper limit on    by 12% and a reduction in the width of the expected one-dimensional confidence interval for   ( 2 ) by 6% (17%), based on a consistent statistical procedure for evaluating the 95% confidence interval as described at the beginning of this section.The observed upper limit on    is reduced by 5%, while the observed one-dimensional confidence interval for   ( 2 ) is increased by 5% (reduced by 16%).
The increase in the width of the observed   confidence interval arises from the fact that this new analysis favours larger, less negative values of the signal strength, corresponding to larger magnitudes of the coupling strength modifier   .The compatibility, considering only statistical uncertainties, between the allowed   interval at 95% CL from this study and that of Ref. [32] is evaluated using a bootstrap technique [100], based on the data events passing the selection of either the previous analysis, or that of the current one, or both.The compatibility between the two results is at the level of 0.3 standard deviations.In each of the regions, a higher category index corresponds to higher BDT scores and more signal-like events.The solid line peaks near 125 GeV are due to single Higgs boson production.Table 3: The expected number of events (estimated by using simulation) from  signals with various   and  2 hypotheses and single Higgs boson production, and the expected number of events from the continuum background, evaluated in the 120 GeV <   < 130 GeV window.For comparison, the number of observed data events is also shown.The uncertainties in the  signals and single Higgs boson backgrounds include the systematic uncertainties discussed in Section 6. Asymmetric uncertainties arise primarily from the theory calculation of the SM ggF  cross-section and the large uncertainty in the yield of single Higgs bosons produced in ggF events in association with heavy-flavour jets, parameterised by a lognormal distribution.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.

Effective field theory interpretation
Anomalous Higgs boson self-interactions or interactions with the other gauge fields and fermions can alter the Higgs boson pair production cross-section and kinematics, as well as the Higgs boson decay rates.The results of the previous section are thus interpreted in the context of two effective field theories to set constraints on the Wilson coefficients of the operators describing these anomalous interactions.
The approach used here follows closely that described in Ref. [34], probing a similar set of operators and benchmark points that directly affect  production: three Wilson coefficients ( ℎℎℎ ,   ℎℎ ,  ℎℎ ) and seven benchmark points [101] of the Higgs effective field theory, and two Wilson coefficients (  ,   ) of the ( † ) 3 and ( † )□( † ) operators of the 'Warsaw' basis [102] of the SM effective field theory.In the SMEFT Lagrangian, the operators O  are multiplied by coefficients   /Λ 2 , where Λ is the energy scale that bounds from above the range of validity of the EFT approach.In this study, a value of Λ = 1 TeV is assumed.In the HEFT interpretation, the only considered coefficient affecting VBF  production is  ℎℎℎ , and thus this production mode is always subdominant relative to ggF .In the SMEFT interpretation, the effects of the operators on VBF  production are similarly expected to be small, since the SMEFT preserves the Higgs doublet structure of the SM and the corresponding cancellation between the   and   diagrams involved in VBF  production.Consequently both interpretations consider only ggF  production while VBF  is assumed to be negligible.

HH bb
Observed 68% CL Observed 95% CL Expected 68% CL Expected 95% CL Best fit SM prediction Figure 6: Likelihood contours at 68% (solid line) and 95% (dashed line) CL in the (  ,  2 ) parameter space, when all other coupling modifiers are fixed to their SM predictions.The corresponding expected contours are shown by the inner and outer shaded regions The SM prediction is indicated by the star, while the best-fit value is denoted by the cross.
Predictions for ggF  production for various values of the Wilson coefficients under study are obtained by applying an event reweighting technique to the SM ggF  sample, similar to the method described in Section 3 to emulate samples with anomalous values of   .The reweighting functions are based on the particle-level    distributions predicted at NLO accuracy in the strong coupling constant for alternative values of the EFT coefficients.For the HEFT interpretation, the functions are taken directly from Ref. [103], while for the SMEFT intepretation, they are computed using Powheg Box v2 with the SMEFT@NLO model [104].For the SMEFT interpretation, a similar reweighting function is also derived for single Higgs boson processes, but instead using the differential distribution of the Higgs boson transverse momentum.
Uncertainties related to PDF,   , and missing higher-order terms in the prediction are included by taking for each analysis category the envelope of the uncertainties from each source, determined with the same procedure as that described in Section 6.In addition, a non-closure uncertainty is estimated by comparing the expected yields from dedicated samples corresponding to specific values of the anomalous couplings to those from the reweighting procedure described above in categories reproducing the analysis selections at generator level.These uncertainties in the expected yield are generally of the order of 10% or less in each category and have a small impact on the results.
In the HEFT interpretation, constraints on the coefficients  ℎℎℎ ,   ℎℎ , and  ℎℎ that describe Higgs boson self-interactions as well as effective  t and  interactions are determined from the data from one-dimensional scans of the profile likelihood function as a function of the coefficients.The operators corresponding to these coefficients do not impact single Higgs boson production and decay at tree level and their effect on the resonant background and on the Higgs boson branching ratios is therefore neglected.
The one-dimensional constraints on the three HEFT coefficients  ℎℎℎ ,   ℎℎ and  ℎℎ are summarised in Table 4.The difference between the  ℎℎℎ constraint and the   constraint previously presented in Figure 5 is mainly due to the lack of VBF production in the former.In addition, the observed constraints are comparable with those of Ref. [34], when evaluated using the same statistical procedure of Ref. [34].The width of the allowed 95% CL interval for  ℎℎ is 20% narrower, while that of the   ℎℎ interval is the same.Figure 7 shows two-dimensional profile log-likelihood contours for the simultaneous variation of the ( ℎℎ ,  ℎℎℎ ) and (  ℎℎ ,  ℎℎℎ ) HEFT coefficients, with the remaining coefficient fixed to its SM value.In addition, upper limits are set on the Higgs boson pair production cross-section for seven benchmark points [101] corresponding to different values of the five coefficients  ℎℎℎ ,   ℎℎ ,  ℎℎ ,  ℎ , and   ℎ , where the latter two correspond to an effective Higgs-gluon interaction and to the Higgs-top Yukawa interaction.The impact of these coefficients on single Higgs boson production and decay is expected to be small compared to the signal and is thus neglected.Defined in Table 5, the benchmark points describe representative signal kinematics and    shape features, and have sensitivities that can vary significantly between one point and another.For example, benchmark 1 results in a very soft    distribution while benchmark 5 produces a more SM-like    distribution with an enhanced tail.
The resulting upper limits on the Higgs boson pair production cross-section through gluon-gluon fusion are shown in Figure 8.For benchmark points 3, 5 and 7, this analysis sets upper limits similar to those set by the search for  → 4 events [34], and, in an analogous way, excludes these scenarios at 95% CL.
The remaining benchmarks (1, 2, 4, and 6) have updated definitions compared to those used in Ref. [34] and therefore the results cannot be directly compared.Benchmark 4 is excluded for the first time at 95% CL by this study, while the other three scenarios are compatible with the data.In the SMEFT interpretation, one-dimensional constraints are derived on the Wilson coefficients after fixing all other coefficients to zero.The results are obtained by including the contributions to the  and  cross-sections from both linear and quadratic terms in the Wilson coefficient expansion.An interpretation in which the expansion is truncated at linear order is poorly constrained due to the dominance of the quadratic term and can yield negative signal cross-sections.The impact of the operators under study on ggF  production parameterised as a function of    and on single Higgs boson production parameterised as a function of the Higgs boson transverse momentum are included in the interpretation.As in the case of HEFT, the coefficients do not impact the Higgs boson decay branching ratios.The one-dimensional constraints on the SMEFT Wilson coefficients in the scenario where the other parameters are fixed to zero, as expected in the SM, are summarised in Table 6.When using the same statistical procedure of Ref. [34] to determine the constraints on the SMEFT Wilson coefficients, the results are only mildly affected, and the size of the 95% CL interval for   (  ) is 38% (10%) smaller than that in Ref. [34].Furthermore, Figure 9 shows two-dimensional likelihood scans as a function of the couplings   and   .

HH bb
Observed limit (95% CL) Expected limit (95% CL) Expected limit ±1 Expected limit ±2 Theory prediction Figure 8: The observed (filled circles) and expected (hollow circles) 95% CL upper limits on the  ggF production cross-section in the SM and for seven HEFT benchmark points defined in Ref. [101].The expected constraints are obtained from a background hypothesis with    = 0.The predicted cross-sections of each of the models under consideration are shown by the crosses.Benchmarks where the filled circles are below the crosses are excluded.The inner and outer shaded bands indicate the ±1 and ±2 variations on the expected limit due to statistical and systematic uncertainties.The contribution from VBF production to the total  production cross-section is neglected.

Conclusion
An updated search for nonresonant Higgs boson pair production in the  b final state is performed using the full Run 2 ATLAS data, corresponding to 140 fb −1 of 13 TeV   collisions.The results supersede and expand upon those of a previous nonresonant search based on the same data sample.Compared to the previous publication, the classification of events in orthogonal event categories is reoptimised to increase the sensitivity to  production in the main production modes, ggF and VBF, and to the Higgs boson self-coupling and quartic coupling to ,  bosons.The sensitivity is increased by 6%-17% depending on the parameter of interest.The statistical procedure for the interpretation of the observed yields in terms of the signal coupling strength modifiers has also been updated.In addition, the results are interpreted in the context of the Higgs and SM effective field theory frameworks to constrain the Wilson coefficients of operators describing anomalous Higgs boson interactions.
No evidence of signal is found.In the most sensitive categories of the analysis a small deficit of events in the signal region leads to a 95% CL upper limit on the  production signal strength    < 4.0 that is lower than the expected value of 5.0 (6.4) in the background-only    = 0 (SM    = 1) hypothesis.
The corresponding observed (expected) one-dimensional intervals at 95% CL for the self-coupling modifier   and the quartic coupling modifier  2 are −1.4 <   < 6.9 (−2.8 <   < 7.8) and −0.5 <  2 < 2.7 (−1.1 <  2 < 3.3), respectively.From these results, one-dimensional limits on the Wilson coefficients of operators affecting Higgs boson pair production in the Higgs effective field theory ( ℎℎℎ ,   ℎℎ ,  ℎℎ ) and SM effective field theory (  ,   ) frameworks are inferred.In the former, the comparison between the predicted gluon-gluon fusion  cross-sections and the corresponding upper limits set by the analysis excludes four of the seven benchmark points considered at 95% CL.While three of these were already excluded by a similar interpretation of the results in the ATLAS search for  production in the 4 final state, one newly proposed benchmark is excluded for the first time by the results presented in this paper.
The one-dimensional constraints on the Wilson coefficients considered in this analysis are up to 38% tighter than those reported previously by ATLAS when evaluated using the same statistical procedure.

Figure 1 :
Figure 1: The Feynman diagrams for the dominant gluon-gluon fusion production processes.In the SM, the (a) trilinear coupling process, (b) box diagram, and the destructive interference between the two processes, contribute to the total cross-section.In the figure,   represents the Higgs boson trilinear coupling modifier.The quark content in the diagram is dominated by the top-quark contribution due to the large top-quark Yukawa coupling to the Higgs boson.The corresponding coupling strength modifier is denoted by   .

Figure 2 :
Figure 2: The VBF production of Higgs boson pairs via (a) the trilinear coupling, (b) the  vertex, and (c) the   production mode.In the figure,   and  2 denote the  and  coupling strength modifiers.

Figure 4 :
Figure 4: Comparison between the diphoton invariant mass distribution in data (points with error bars) and the background-only fit (solid line) for the four low mass (a-d) and three high mass (e-g) categories of the  →  b search.In each of the regions, a higher category index corresponds to higher BDT scores and more signal-like events.The solid line peaks near 125 GeV are due to single Higgs boson production.

Figure 5 :
Figure 5: Observed (solid line) and expected (dashed line) value of −2 ln Λ as a function of (a)   and (b)  2 , when all other coupling modifiers (including, respectively,  2 or   ) are fixed to their SM predictions.

Figure 7 :
Figure7: Likelihood contours at 68% (solid line) and 95% (dashed line) CL in the (a)  ℎℎ versus  ℎℎℎ and (b)   ℎℎ versus  ℎℎℎ HEFT parameter space, with the remaining coefficient fixed to its SM value.The corresponding expected contours are shown by the inner and outer shaded regions.The SM prediction is indicated by the star, while the best-fit value is denoted by the cross.

Figure 9 :
Figure9: Likelihood contours at 68% (solid line) and 95% (dashed line) CL in the   versus   SMEFT parameter space.The corresponding expected contours are shown by the inner and outer shaded regions.The SM prediction is indicated by the star, while the best-fit value is denoted by the cross.
Higgs boson decaying into  b and the other one into .In addition to the samples in Table1, a ggF  sample was generated with the same settings as the nominal sample but with the non-SM value of the self-coupling modifier   = 10, and then passed through the detector simulation and the reconstruction algorithms.A reweighting technique based on the particle-level invariant mass    of the Higgs boson pair is applied to the   = 1 sample to determine the ggF  signal yield and kinematic distributions for any value of Higgs boson pair production, especially    , are significantly affected by the values of   and  2 .In particular, ggF and VBF  production with values of   close to the SM expectation lead to rather large values of    , while for   significantly different from one the  invariant mass spectrum is relatively soft.Anomalous values of  2 also lead to events, produced via VBF, with a large invariant mass of the Higgs boson pair.The events are therefore classified in two regions based on the modified four-body invariant mass  *  b =   b − (  b − 125 GeV) − (  − 125 GeV): a high mass ( *  b > 350 GeV) region and a low mass ( *  b ≤ 350 GeV) region.The use of  *  b over   b

Table 4 :
The observed and expected 95% CL constraints on the HEFT Wilson coefficients, obtained from onedimensional scans of the profile log-likelihood assuming that all other Wilson coefficients are fixed to their SM values.The contribution from VBF  production is subdominant to that from ggF and is neglected.

Table 6 :
The observed and expected 95% CL constraints on the SMEFT Wilson coefficients, obtained from one dimensional scans of the profile log-likelihood assuming that all other Wilson coefficients are fixed to their SM values.The contribution from VBF production is neglected.