Search for Higgs boson pair production in the $\gamma \gamma WW^*$ channel using $pp$ collision data recorded at $\sqrt{s}$ $=13$ TeV with the ATLAS detector

Abstract

Searches for non-resonant and resonant Higgs boson pair production are performed in the $\gamma \gamma WW^*$ channel with the final state of $\gamma \gamma \ell \nu jj$ using 36.1 $\hbox {fb}^{-1}$ of proton–proton collision data recorded at a centre-of-mass energy of $\sqrt{s}=$13 TeV by the ATLAS detector at the Large Hadron Collider. No significant deviation from the Standard Model prediction is observed. A 95% confidence-level observed upper limit of 7.7 pb is set on the cross section for non-resonant production, while the expected limit is 5.4 pb. A search for a narrow-width resonance X decaying to a pair of Standard Model Higgs bosons HH is performed with the same set of data, and the observed upper limits on $\sigma (pp\rightarrow X)\times B(X\rightarrow HH)$ range between 40.0 and 6.1 pb for masses of the resonance between 260 and 500 GeV, while the expected limits range between 17.6 and 4.4 pb. When deriving the limits above, the Standard Model branching ratios of the $H{\rightarrow }{\gamma \gamma }$ and $H{\rightarrow }WW^*$ are assumed.

1 Introduction

A particle consistent with the Standard Model (SM) Higgs boson (H) was discovered by both the ATLAS and CMS experiments at the Large Hadron Collider (LHC) in 2012 [1, 2]. Various studies of its properties have been performed [3,4,5,6,7], and no significant deviation from the SM predictions has been found. The SM Higgs boson is a strong probe of physics beyond the SM. This paper documents searches for both non-resonant and resonant production of Higgs boson pairs ($HH$) in the semileptonic $\gamma \gamma WW^*$ final state using 36.1 $\hbox {fb}^{-1}$ of proton–proton ($pp$) collision data recorded by the ATLAS detector at a centre-of-mass energy of $\sqrt{s} = 13\,\hbox {TeV}$. Previous searches for Higgs boson pair production have been performed by both the ATLAS and CMS experiments with data recorded at $\sqrt{s} = 8\,\hbox {TeV}$ in the final states $b\bar{b}b\bar{b}$  [8], $b\bar{b}\gamma \gamma$  [9, 10], $b\bar{b}\tau ^{+}\tau ^{-}$  [11,12,13] and $\gamma \gamma WW^*$  [11], as well as multi-lepton and multi-photon [14]. The $pp$ collision data at $\sqrt{s} = 13\,\hbox {TeV}$ have been analysed in order to search for Higgs boson pairs in the final states $b\bar{b}b\bar{b}$  [15] and $b\bar{b}WW^*$  [16]. No significant excess was observed compared to the SM prediction. However, it is important to explore the $13\,\hbox {TeV}$ data in the channels that are not covered yet, such as the $\gamma \gamma WW^*$ channel presented in this paper. Although this decay channel is not the most sensitive amongst all possible Higgs boson decays, it relies on the Higgs boson couplings to vector bosons, which are already relatively well measured. Furthermore, this channel will contribute to the final combination of all measurable HH decays.

The SM prediction of the Higgs boson pair production cross section is several orders of magnitude smaller than the single-Higgs-boson production rate [17], due to additional ttH or HHH vertices, an additional on-shell Higgs boson that reduces the kinematic phase space, and the fact that the two leading-order (LO) Feynman diagrams have strong destructive interference [18]. In Fig. 1a, the so-called box diagram represents Higgs boson pair production via a heavy-quark loop, where the cross section scales with the squared value of the ttH or bbH coupling constants. In Fig. 1b, the so-called triangle diagram contributes to Higgs boson pair production via the exchange of a virtual Higgs boson and is the only tree-level diagram sensitive to the Higgs boson self-coupling constant ($\lambda _{HHH}$), the squared value of which scales the cross section.

In many beyond-the-SM (BSM) scenarios, Higgs boson pair production can be enhanced by modifying the ttH, bbH or $\lambda _{HHH}$ coupling constants from their SM values, reducing the effect of the destructive interference [19] between Fig. 1a and Fig. 1b, or by replacing the virtual Higgs boson with an intermediate scalar resonance, cf. Fig. 1c. Various BSM models with extended Higgs sectors predict a heavy Higgs boson decaying into a pair of Higgs bosons similar to the one in the SM. Such models include the two-Higgs-doublet models (2HDM) [20], the minimal supersymmetric extension of the SM [21], twin Higgs models [22] and composite Higgs models [23, 24]. Heavy resonances, other than heavy Higgs bosons, that can decay into a pair of SM Higgs bosons, are predicted in different models, and could for instance be gravitons [25], radions [26] or stoponium [27].

This paper reports searches for non-resonant and resonant production of pairs of Higgs bosons in the semileptonic $\gamma \gamma WW^*$ final state ($\gamma \gamma \ell \nu jj$), i.e. with two photons, two jets, one charged lepton and a neutrino. This final state benefits from the large branching fraction of $H{\rightarrow }WW^*$  [17], a characteristic signature from two photons and one lepton, as well as the excellent resolution of the diphoton invariant mass $m_{\gamma \gamma }$, which provides good discrimination from a smooth continuum background composed of multi-photon and multi-jet SM processes. Given the expected sensitivity in $13\,\hbox {TeV}$ data, the di-Higgs-boson mass range between 260 and 500 GeV is explored in the search for a scalar resonant Higgs boson pair production.

Fig. 1

Feynman diagrams for leading-order Higgs boson pair production in the SM through a a heavy-quark loop, b the Higgs self-coupling, and c an intermediate heavy resonance in a BSM scenario. The total SM contribution is the sum of the two modes depicted in a and b, which have significant destructive interference. Physics beyond the SM can enhance Higgs boson pair production either by modifying the Higgs boson coupling constants from their SM values in a and/or b, or by an additional s-channel exchange of an intermediate scalar resonance in c

2 The ATLAS detector

The ATLAS experiment [28] is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and nearly $4\pi$ coverage in solid angle.¹ It consists of an inner tracking detector (ID) surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic (EM) and hadronic calorimeters, and a muon spectrometer (MS). The ID covers the pseudorapidity range $|\eta | < 2.5$ and consists of silicon pixel, silicon microstrip, and transition-radiation tracking systems. The innermost pixel layer, the insertable B-layer [29], was installed at a mean radius of 3.3 cm after Run 1, and has been operational since the beginning of Run 2. Lead/liquid-argon (LAr) EM sampling calorimeters with high granularity provide energy measurements of EM showers. A steel/scintillator-tile hadronic calorimeter covers the central pseudorapidity range ($|\eta | < 1.7$). The endcap and forward regions are covered by LAr calorimeters for EM and hadronic energy measurements up to $|\eta | = 4.9$. The MS surrounds the calorimeters and is based on three large air-core toroid superconducting magnets with eight coils each and with bending power in the range 2.0–7.5 T m. It includes a system of fast detectors for triggering purposes and precision tracking chambers. A dedicated two-level trigger system is used to select events [30]. The first-level trigger is implemented in hardware and uses a subset of the detector information to reduce the accepted event rate to at most 100 kHz. This is followed by a software-based trigger level that reduces the accepted event rate to an average of 1 kHz.

3 Data and simulated samples

3.1 Data samples

The full set of $pp$ collision data collected during 2015 and 2016 are used in this analysis. The two datasets were recorded at the same centre-of-mass energy $\sqrt{s}$ $=$13 TeV, albeit with different beam conditions. Beam intensities in 2016 were typically higher than in 2015, resulting in a higher instantaneous luminosity and a larger number of $pp$ collisions in each bunch crossing. The integrated luminosity of the combined 2015+2016 dataset used in this analysis is $36.1\pm 0.8\,\hbox {fb}^{-1}$. This dataset were collected in run periods during which all subsystems were operational. The events are collected with a trigger requiring the presence of at least two photons, one with a transverse energy $E_{\text {T}}$ > 35 GeV and the second with $E_{\text {T}}$ > 25 GeV, and the longitudinal and transverse profiles of the EM shower were required to be consistent with those expected for a photon. The corresponding trigger efficiency reaches about 99% for the events that pass the event selection of the analysis.

3.2 Simulated event samples

Simulated Monte Carlo (MC) samples are used to estimate the signal acceptance and study the modelling for both non-resonant SM Higgs boson pair production and resonant BSM Higgs boson pair production. MC samples are also used to estimate the acceptance and study the modelling for SM single-Higgs-boson production processes, and to study the modelling of the SM continuum background from events with multiple photons and jets (Sect. 5), which is the dominant background in the analysis. Eventually, it is estimated by a data-driven method for both its normalisation and shape.

Table 1 Simulated signal samples

The simulated samples for signals are listed in Table 1. The event generator MadGraph5_aMC@NLO 2.2.3 [31] was used for the production of non-resonant [32] and resonant [33] signal MC samples at next-to-leading order (NLO) in QCD, where four values of the resonance mass ($m_{X}$ = 260, 300, 400 and 500 GeV) are considered. The events were generated by a Higgs Effective Field Theory (HEFT) using the MC@NLO method [34] and were reweighted in order to take into account the effects of the finite top-quark mass. The parton shower was implemented using Herwig++ 2.7.1 [35] with a set of tuned underlying-event parameters called the UEEE5 tune² [36], and the parton distribution function (PDF) set CTEQ6L1 [37] was used.

• For non-resonant Higgs boson pair production, the inclusive cross section s are normalised to the SM prediction of 33.41 fb [17, 38], calculated at NNLO in QCD, including resummation of soft-gluon emission at next-to-next-to-leading-logarithmic (NNLL) accuracy, as prescribed by the LHC Higgs Cross Section Working Group [17]. The effect of the finite top-quark mass is also taken into account at NLO [39].

• For resonant Higgs boson pair production, a narrow decay width, which is negligible compared to the experimental mass resolution, is assumed. The interference between non-resonant and resonant Higgs boson pair production is implemented in the generator. The interference is minimal and remains negligible when a narrow decay width is assumed.

Table 2 Simulated SM single-Higgs-boson background samples with $m_{H} =$ 125 GeV

Table 2 lists the simulated samples for the dominant SM single-Higgs-boson production modes: gluon–gluon fusion (ggF), vector-boson fusion (VBF), associated production with a W or Z boson ($VH$), and associated production with a pair of top quarks ($t\bar{t}H$). For all these processes, the Pythia 8.186 parton shower is used for the modelling of non-perturbative effects. The AZNLO tune [40] is used in ggF, VBF and $VH$ simulations, while the A14 tune is used in $t\bar{t}H$ simulations.

• gluon–gluon fusion:  The ggF production is accurate to NNLO in QCD, using the Powheg method [41] for matching the matrix element with the parton shower, and the MiNLO method [42, 43] to simultaneously achieve NLO accuracy for inclusive Higgs boson production. Furthermore, a reweighting procedure was performed using the HNNLO program [44,45,46] to achieve full NNLO accuracy [47]. This sample is referred to as NNLOPS. The PDF4LHC15 NLO PDF set [48] was used. The inclusive cross section of the ggF production is normalised to the calculation at next-to-next-to-next-to-leading-order (N$^3$LO) QCD and NLO electroweak (EW) accuracies [49].

• VBF and VH:  VBF and $VH$ production was simulated at NLO in QCD with Powheg-Box v2 [41, 50, 51] using the PDF4LHC15 NLO PDF set. The inclusive VBF contribution is normalised to the cross section calculated with NLO QCD and NLO EW corrections [52,53,54] with an approximate NNLO QCD correction applied [55]. The contributions are normalised to cross section s calculated with NNLO QCD [56] and NLO EW corrections [57] for WH and $q\bar{q}\rightarrow ZH$ and at NLO and next-to-leading-logarithm (NLL) accuracy in QCD for $gg\rightarrow ZH$ [58].

• ${\varvec{t}}\bar{{\varvec{t}}}{\varvec{H}}$:  The $t\bar{t}H$ production is simulated using MadGraph5_aMC@NLO 2.2.3 and its inclusive cross section is normalised to a calculation with NLO QCD and NLO EW corrections [59,60,61,62].

Processes of continuum backgrounds of multiple photons and jets with either one or zero leptons were simulated with MadGraph5_aMC@NLO 2.2.2, interfaced with the parton shower model in Pythia 8.186.

Multiple $pp$ collisions in each bunch crossing, “pile-up”, were simulated with the soft QCD processes of Pythia 8.186 using the A2 tune [63] and the MSTW2008LO PDF set [64]. An additional event-level reweighting is performed in order to ensure that the distribution of the average number of interactions per bunch crossing matches that occurring in the data used in this analysis. The particles in the final states of the generated processes were passed through either a Geant4  [65] simulation of the ATLAS detector, or through the ATLAS fast simulation framework [66], which has been extensively validated against the Geant4 simulation model. The output from the detector simulation is then analysed using the same reconstruction software as the data. The MC samples for single-Higgs-boson production were simulated with the Geant4 framework, while the other samples used in this analysis were produced with the ATLAS fast simulation framework.

4 Object and event selection

The event selection is based on the properties of the visible objects in the final state, which includes one charged lepton (electron or muon), two jets, and two photons. These objects are reconstructed from detector-level objects, such as energy clusters in the EM calorimeter and tracks in the ID, as well as in the MS.

4.1 Object reconstruction

Photon candidates are reconstructed from clusters of energy deposited in the EM calorimeter [67]. If the candidates are matched with a reconstructed conversion vertex or tracks consistent with the hypothesis of a $\gamma {\rightarrow }e^{+}e^{-}$ conversion, they are classified as converted photon candidates.³ If the matched track is consistent with the hypothesis of an electron produced in the beam interaction region, they are classified as electron candidates. If the candidates are not matched with a reconstructed conversion vertex or tracks satisfying the conversion requirement, they are classified as unconverted photon candidates. The energy is determined by summing the energies of all cells that belong to the associated cluster [68] and is corrected using a combination of simulation-based and data-driven calibration factors [69] determined from $Z\rightarrow e^+e^-$ events collected in 2015 and 2016. The photon energy resolution in simulation is corrected to match the resolution in data [67]. The reconstructed photon candidates are required to meet “tight” photon identification criteria [68], which are based on the lateral and longitudinal energy profiles of EM showers in the calorimeter. The identification efficiency is measured as a function of the transverse energy of photons ($E_{\text {T}}^{\gamma }$). It ranges from 90 to 98% for converted photons and from 85 to 95% for unconverted photons, in the $E_{\text {T}}^{\gamma }$ interval between 25 and 200 $\text {GeV}$. To suppress the background from jets misidentified as photons, all reconstructed photon candidates are required to meet a set of calorimeter- and track-based isolation criteria [70]. A calorimeter-based isolation variable $E^{\text {iso}}_{\text {T}}$ is defined as the sum of the transverse energies ($E_{\text {T}}$) of all positive-energy topological clusters of calorimeter cells [71] within $\Delta R = 0.2$ of the photon candidate, excluding the energy of the photon candidate itself. The selection applied to the calorimeter-based isolation variable is $E^{\text {iso}}_{\text {T}} < 0.065E_{\text {T}} ^{\gamma }$. A track-based isolation variable $p^{\text {iso}}_{\text {T}}$ is defined as the scalar sum of the transverse momenta ($p_{\text {T}}$) of tracks with $p_{\text {T}}$ > 1 GeV within $\Delta R = 0.2$ of the photon candidate, excluding tracks from photon conversions. The selection on the track-based isolation variable is $p^{\text {iso}}_{\text {T}} < 0.05E_{\text {T}}^{\gamma }$. Only photon candidates with $|\eta |$ < 2.37 are considered, excluding the transition region between the barrel and endcap calorimeters (1.37 < $|\eta |$ < 1.52).

Electron candidates are reconstructed from clusters of energy deposited in the EM calorimeter matched to a track in the inner detector, as described above. A likelihood-based (LH) algorithm is used [72] to perform the electron identification against the background from jets or non-prompt electrons. Electron candidates are identified according to the “medium LH” criteria. Muon candidates are identified by matching a reconstructed ID track with a reconstructed MS track [73]. The identification classifies muon candidates as either “loose” or “medium”, based on the number of hits in the different ID and MS subsystems, and on the significance of the difference $|q/p_{\text {MS}}-q/p_{\text {ID}}|$, where q is the charge and p is the momentum of the muon candidate, as well as on the energy deposit in the tile hadronic calorimeters. The “medium” candidates are used in the analysis. An efficiency ranging from 84 to 93% as a function of $E_{\text {T}}$ or $p_{\text {T}}$ is achieved in the combined identification and reconstruction of electrons, and 96% (above 98%) in muon identification (reconstruction), in the range where the objects are selected. The electron (muon) is required to pass the “Loose” (“GradientLoose”) isolation criterion based on the sum of $p_{\text {T}}$ of tracks lying within a cone of $\Delta R=$ min (10 $\hbox {GeV/}p_{\text {T}}^{e (\mu )}$, 0.3) and the sum of $E_{\text {T}}$ of topological clusters of calorimeter cells within a cone of $\Delta R=$ 0.2 (0.2) around the electron (muon) candidate, excluding the contributions from the electron (muon) candidate. With these requirements the isolation efficiencies for electrons (muons) are above 99% ($0.057p_{\text {T}}^{\mu }+95.57\%$) [74, 74]. Finally, the electron candidates are required to have $E_{\text {T}} > 10$ GeV and $|\eta |<2.47$, excluding the transition region between the barrel and endcap calorimeters (1.37 < $|\eta |$ < 1.52), whereas the muon candidates are required to have $p_{\text {T}} > 10$ GeV and $|\eta |<2.7$.

Jets are reconstructed via the FastJet package [74] using the anti-$k_{t}$ clustering algorithm [75] with a radius parameter $R=$ 0.4. The jet energies are determined at the EM scale and calibrated using particle-level correction factors based on a combination of simulation and data [76,77,78,79,80]. Jets are required to have $|\eta |$ < 2.5 and $p_{\text {T}}$ > 25 GeV. In addition, a jet-vertex tagging algorithm (JVT)  [81] is applied to jets with $|\eta |$ < 2.4 and $p_{\text {T}}$ < 60 GeV in order to suppress jets originating from pile-up interactions. In this algorithm, a multivariate discriminant based on two track-based variables is constructed to reject pile-up jets while maintaining a high efficiency for the hard-scatter jet independent of the number of primary vertices in the event. The selected jets are classified as $b\text {-jet}$ s using a multivariate technique [82, 83], which takes advantage of the information about secondary vertices, the impact parameters of the associated tracks and the topologies of decays of heavy-flavour hadrons. The b-tagging working point is selected to have an efficiency of 70% for a b-jet from $t\bar{t}$ decays, with a rejection factor of 12 for jets originating from c-quarks (c-jets), and of close to 400 for jets initiated by light-flavour quarks or gluons (light-flavour jets).

An overlap removal procedure is performed in the following order to avoid double counting of detector-level objects when reconstructing physics objects. Electrons with $\Delta R(e,\gamma )<$ 0.4 are removed. Jets with $\Delta R(\text {jet},\gamma )<$ 0.4 or $\Delta R(\text {jet},e)<$ 0.2 are removed. Electrons with $\Delta R(e,\text {jet})<$ 0.4 are removed. Muons with $\Delta R(\mu ,\gamma )<$ 0.4 or $\Delta R(\mu ,\text {jet})<$ 0.4 are removed.

4.2 Event selection

The events passing the diphoton trigger are required to contain at least two jets, no $b\text {-jet}$, and at least one charged lepton (e or $\mu$, but including contributions from fully leptonic $\tau$-lepton decays) in the final state. The two photon candidates with the leading (sub-leading) $E_{\text {T}}$ are required to satisfy $E_{\text {T}} ^{\gamma }/m_{\gamma \gamma }>$ 0.35 (0.25). The $b\text {-jet}$ veto suppresses the $t\bar{t}H$ process. Furthermore, the transverse momentum of the diphoton system ($p_{\text {T}} ^{\gamma \gamma }$) is required to be larger than 100 GeV for maximising the sensitivity and keeping at least 70% of signal events. This requirement suppresses continuum background events when searching for non-resonant Higgs boson pair production, or resonant production with resonance masses of 400 GeV or higher. However, the $p_{\text {T}} ^{\gamma \gamma }$selection is omitted in the search for resonance masses below 400 GeV due to a limited separation between signal and continuum background in this kinematical region, as can be seen in Fig. 2. These final selection criteria, together with a requirement on the invariant diphoton mass of $105\,\hbox {GeV}< m_{\gamma \gamma } < 160\,\hbox {GeV}$, define the event sample on which the signal search is performed for the various assumed signal models. A data “sideband” sample is selected applying the same criteria, but excluding the Higgs mass region $m_{\gamma \gamma }$ 121.7–128.5 GeV, and can be used together with other samples to study the continuum background.

If there were an observable signal, one of the Higgs bosons would be directly visible in the $m_{\gamma \gamma }$ distribution. The combination of two jets and at least one charged lepton would be consistent with $H{\rightarrow }WW^*$ for the other Higgs boson. Its signature would in principle be enhanced by a missing transverse energy ($E_{\text {T}}^{\text {miss}}$) requirement to indicate a neutrino, but a selection on $E_{\text {T}}^{\text {miss}}$ was found not to produce any significant improvement in sensitivity, and so was not applied. The magnitude of $E_{\text {T}}^{\text {miss}}$  [84, 85] is measured from the negative vectorial sum of the transverse momenta of all photon, electron and muon candidates and of all hadronic jets after accounting for overlaps between jets, photons, electrons, and muons, as well as an estimate of soft contributions based on tracks.

Fig. 2

Distributions of the reconstructed transverse momenta of the diphoton system with all event selections, except the $p_{\text {T}} ^{\gamma \gamma }$selection, applied for various signal models, as well as sideband data, normalised to unit area

After all selections described above, the combined acceptance and selection efficiency for non-resonant production is 8.5%, while it ranges from 6.1 to 10% as a function of the mass of the resonance ($m_X$) from 260 to 500 GeV, as shown in Table 3. The efficiency for the non-resonant Higgs boson pair production is at the same level as the efficiency for the high-mass resonant production, as the Higgs bosons and their decay products tend to exhibit large transverse momenta due to the box diagram shown in Fig. 1a.

Table 3 The combined acceptance and efficiency for non-resonant and resonant with different scalar resonance masses $m_{X}$, with and without a $p_{\text {T}} ^{\gamma \gamma }$selection

5 Signal and background estimation

A fit to the $m_{\gamma \gamma }$ distribution is performed to extract the signal yield as described in Sect. 7. The shapes of both the signal and background distributions are modelled with analytical functions. For both Higgs boson pair production and single-Higgs-boson processes, the $m_{\gamma \gamma }$ distributions are modelled with double-sided Crystal Ball functions [86]. Their shape parameters are determined by a fit to simulated samples. The single-Higgs-boson contribution is normalised to the SM cross-sections as described in Sect. 3.2. Higgs boson pair production is regarded as a background to the resonant search. Its contribution is also set to the SM prediction of Sect. 3.2.

The continuum background is modelled with an exponential function of a second-order polynomial. Several functional forms were evaluated by fitting the sidebands in data and MC samples under different conditions of photon purity and lepton multiplicity. Photon purity was lowered, compared to the final data selection, by reversing the requirements on photon isolation or identification. For higher photon purity, MC samples with prompt photons were used. The lepton multiplicity was varied to be zero or at least one. For all combinations of conditions, the exponential function with a second-order polynomial gave the best fits, with satisfactory $\chi ^2$, and was chosen to model the continuum background. The shape parameters and normalisation are free to float in the final fit to the data. Since any functional form might introduce spurious signals, this effect is estimated with a sample mixing irreducible prompt-photon background from simulation and reducible backgrounds from data, as described in Sect. 6.

The expected numbers of signal and background events are shown in Table 4 together with the number of events observed in data. Only events within a mass window of $m_{H} \pm 2$ $\sigma _{m_{\gamma \gamma }}$ are reported, where the Higgs boson mass ($m_{H}$) is taken to be 125.09 GeV [87] and the diphoton mass resolution ($\sigma _{m_{\gamma \gamma }}$) is 1.7 GeV and is obtained from simulation. The dominant background is from continuum processes with multiple photons and jets. A small background arises from SM single-Higgs-boson production processes, among which $t\bar{t}H$ and WH productions give the leading contributions with, respectively, a fraction of 41.5% (39.2%) and 23.3% (22.5%) of the whole single-Higgs-boson contribution with (without) the $p_{\text {T}} ^{\gamma \gamma }$> 100 GeV selection.

Table 4 Numbers of expected and observed events in the $m_{H} \pm 2$ $\sigma _{m_{\gamma \gamma }}$ mass window with or without a $p_{\text {T}} ^{\gamma \gamma }$selection. A cross section of 33.41 fb is assumed for non-resonant Higgs boson pair production when it is considered as a background in resonant searches. The resulting yields are determined from the fit to data by integrating the resulting functional forms over the selected $m_{\gamma \gamma }$ range. The error in each yield includes both the statistical and systematic uncertainties, as discussed in Sect. 6

6 Systematic uncertainties

6.1 Theoretical uncertainties

Theoretical uncertainties in the prediction of the cross section of single Higgs bosons are estimated from variations of the normalisation and factorisation scales, PDF, and the running QCD coupling constant ($\alpha _\text {S}$) [17]. Among the dominant production modes $t\bar{t}H$ and $VH$, the cross section of $t\bar{t}H$ has the largest uncertainty: up to 9.2% in the scale variations, up to 3.0% in the PDF variations, and 2.0% in the $\alpha _\text {S}$ variations, as prescribed by the LHC Higgs Cross Section Working Group [17].

The theoretical uncertainties in the efficiency times acceptance ($\epsilon \times A$) are estimated from scale, PDF and parton shower variations. The scale uncertainty ranges from 2.1 to 4.1% for resonant production and is 3.4% for non-resonant production. The PDF uncertainty is around 2.5 and 3.0% for the resonant and non-resonant production, respectively. The parton shower uncertainty is estimated by comparing Pythia 8 and Herwig++ as two different shower models, and ranges from 6.0% at $m_{X} = 500$ GeV to 29.6% at $m_{X} = 260$ GeV for resonant production, and is 7.8% for non-resonant production. This uncertainty is large in low-mass resonant production because the jet spectrum at low-$p_\text {T}$ is more susceptible to variations in the parton shower model. Non-resonant Higgs boson pair production is considered as a background in the search for resonant Higgs boson pair production. The scale, PDF, $\alpha _\text {S}$ and HEFT uncertainties in the calculation of the cross section for SM Higgs boson pair production are also taken into account. These values are 6.0%, 2.1%, 2.3%, and 5.0%, respectively, following the recommendations in Ref. [17]. Further uncertainties arising from the $H{\rightarrow }{\gamma \gamma }$ and $H{\rightarrow }WW^*$ branching ratios (B) are considered as well. They are 2.1% and 1.5% [17], respectively.

6.2 Modelling uncertainties in the continuum background

The exponential function of a second-order polynomial is determined to provide the simplest and most robust functional form for modelling the continuum background as described in Sect. 5. The uncertainties in the modelling are estimated by fitting a signal-plus-background model to a simulated background-only sample that has such a large number of events that its own statistical uncertainty does not affect the test results. The fitted number of signal events ($n_\text {ss}$) quantifies spurious signal events. The fits are performed with the assumed $m_H$ ranging from 120 to 130 GeV in steps of 0.5 GeV. The maximum value of the fitted signal yields |$n_\text {ss}$| is regarded as a bias in the yields due to the background modelling (the spurious signal), and is, conservatively, taken into account in the fit as the modelling uncertainty. The fitted |$n_\text {ss}$| value reaches as large as 0.46 when not applying the $p_{\text {T}} ^{\gamma \gamma }$ selection, and 0.26 when applying the selection. The simulated background-only samples include the irreducible process of $\gamma \gamma \ell \nu jj$ and the reducible processes represented by events where one or two hadronic jets are misidentified as photons. The reducible processes are modelled by the data events with reversed photon identification or isolation requirements. The two components are combined according to the measured diphoton purity, which is about 88% (90% with $p_{\text {T}} ^{\gamma \gamma }$selection) and normalised according to the number of selected data events.

6.3 Experimental uncertainties

The uncertainty in the measurement of the combined 2015+2016 integrated luminosity is 2.1%. It is derived, following a methodology similar to that detailed in Ref. [88], from a calibration of the luminosity scale using x–y beam-separation scans performed in August 2015 and May 2016. All processes that are estimated using simulation are affected by the uncertainty in the luminosity measurement.

The efficiency of the diphoton trigger is estimated using bootstrap methods [89] with a systematic uncertainty of 0.4%. The photon identification uncertainty is obtained by varying the data-to-simulation efficiency corrections within their uncertainties, derived from control samples of photons from radiative Z boson decays and from inclusive $\gamma$ events, and of electrons from $Z\rightarrow e^+e^-$ decays. A maximal uncertainty of 1.7% in the yields is evaluated in all of the SM single-Higgs-boson, SM di-Higgs-boson and BSM Higgs boson production processes. The photon–track isolation uncertainty is derived from measurements of the uncertainty in the data-to-simulation efficiency corrections using inclusive-photon control samples, while the uncertainty from the calorimeter isolation requirement is evaluated from the difference between applying and not applying corrections derived from inclusive-photon events to the calorimeter isolation variable in the simulation. In general, the overall isolation uncertainty is less than 1%. The uncertainties from the photon energy resolution and scale affect the yields by less than 0.2%. The relevant impact on the shape of the diphoton invariant mass is also considered by introducing variations of the resolution and mean values of the fit function and is estimated using simulation. The photon energy resolution varies the resolution of the $m_{\gamma \gamma }$ shape by 5.2–11.4%, while the photon energy scale affects the mean value by about 0.5%. The jet energy scale (JES) and the corresponding uncertainties are derived from both simulation and in situ calibration using data [77, 90]. This affects the event selection efficiency by 2.4–9.9%, depending on the process. The jet energy resolution (JER) uncertainty is evaluated by smearing jet energies according to the systematic uncertainties of the resolution measurement [80, 91], and its impact on the event selection efficiency ranges from 0.1 to 1.6%. The b-tagging uncertainties is derived separately for b-jets, c-jets and light-flavour jets [82]. Overall, their impact on the yields is not more than 4%. Uncertainties arising from the reconstruction, identification and isolation of both the electron and muon candidates [72, 73], are propagated to the event yield variations, and they are found to have an impact of less than 1%. Finally, the pile-up reweighting procedure, which matches the distribution of the number of interactions per bunch crossing between simulation and data, has associated systematic uncertainties of less than 1%. All experimental uncertainties are correlated among all processes that use simulation to model the yields and the kinematics. A summary of the systematic uncertainties in the expected yields of the di-Higgs-boson and single-Higgs-boson production is presented in Table 5. In the search for non-resonant Higgs boson pair production, SM Higgs boson pair production is considered to be the signal process, while single Higgs boson production is considered to be a background. In the search for the resonant Higgs boson pair production, both SM single Higgs boson production and non-resonant Higgs boson pair productions are considered to be background processes.

Table 5 Summary of relative systematic uncertainties, in percent, propagated to the yields for the MC-estimated processes. Entries marked by ‘–’ indicate that the systematic uncertainty is not applicable for the corresponding process. The extrapolation uncertainties in b-tagging include two components: one is from the extrapolation to high-$p_{\text {T}}$ ($p_{\text {T}}$ $>300$ GeV) jets and the other one is from extrapolating c-jets to $\tau$-jets. The values for resonant production shown here assume $m_X$ = 260 GeV. Several theoretical uncertainties are reported for the cross section ($\sigma$) and the combined efficiency and acceptance ($\epsilon \times A$)

7 Results

A fit to $m_{\gamma \gamma }$ is performed in the signal region to extract the signal yield. The statistical model is constructed with a likelihood function:

$$\begin{aligned} \mathcal {L}(\mu ,\varvec{\theta })&= \prod _{i} \Big ( ( n_{\mathrm {Signal}}(\mu ,\varvec{\theta }) + n_{\mathrm {ss}} ) \times f^{1}_{\mathrm {DSCB}}\left( m_{\gamma \gamma }^i,\varvec{\theta } \right) + n_{\mathrm {Cont}} \nonumber \\&\quad \times f_{\mathrm {Cont}}\left( m_{\gamma \gamma }^i,\varvec{\theta } \right) + n_{\text {SM-one-Higgs}}(\varvec{\theta }) \times f^{2}_{\mathrm {DSCB}}\left( m_{\gamma \gamma }^i,\varvec{\theta } \right) \nonumber \\&\quad + n_{\text {SM-di-Higgs}} \times f^{3}_{\mathrm {DSCB}}\left( m_{\gamma \gamma }^i,\varvec{\theta } \right) \Big ) \prod G(0| \varvec{\theta} 0, 1) \end{aligned}$$

(1)

• i stands for the event index,

• $n_{\mathrm {Signal}}$ is the expected number of signal events,

• $\mu$ is the cross section (times the branching fraction of $X\rightarrow HH$) of non-resonant (resonant) production,

• $n_{\mathrm {ss}}$ is the estimated spurious signal yield due to our choice of continuum background modelling,

• $f^{1}_{\mathrm {DSCB}}$ is the probability density function (pdf) of a double-sided Crystal Ball distribution for signal,

• $n_{\mathrm {Cont}}$ is the expected number of continuum background events,

• $f_{\mathrm {Cont}}$ is the pdf of the continuum background, i.e. an exponential function of a second-order polynomial,

• $n_{\text {SM-one-Higgs}}$ is the expected number of single-Higgs-boson events, which is set to the SM prediction and can vary with uncertainties,

• $f^{2}_{\mathrm {DSCB}}$ is the pdf of a double-sided Crystal Ball distribution for the SM single-Higgs-boson background,

• $n_{\text {SM-di-Higgs}}$ is the expected number of the SM di-Higgs-boson events,

• $f^{3}_{\mathrm {DSCB}}$ is the pdf of a double-sided Crystal Ball distribution for SM di-Higgs-boson background,

• $G(0|\varvec{\theta } 1)$ is the pdf of a Gaussian distribution used to constrain the nuisance parameters $\varvec{\theta }$ that model systematic uncertainties as introduced in Sect. 6.

Equation (1) is used directly for the BSM resonant signal searches. For the non-resonant SM Higgs boson pair search, the SM Higgs boson pair term is removed.

Fig. 3

Invariant mass spectrum of the diphoton system in the searches for both resonant and non-resonant Higgs boson pair production, with the corresponding backgrounds for a $m_{X} = 260$ GeV without any $p_{\text {T}} ^{\gamma \gamma }$selection and b the non-resonant case with a $p_{\text {T}} ^{\gamma \gamma }$> 100 GeV selection. Fits to $m_{\gamma \gamma }$ are performed using the full signal-plus-background model. In each plot, only the background component is present. The shape parameters and normalisation of the continuum background model are determined in the fits. The “SM Higgs boson” in a contains the single-Higgs-boson background and SM di-Higgs-boson background. The band shows the uncertainty of the “Total background” in the upper panel and is calculated by a sampling method. The bottom panel shows the difference between the number of events in data and the estimated number of background events, as determined by the fits

The distributions in the final signal-plus-background fit using the likelihood function above are shown for two sets of selections separately: in Fig. 3a without requiring the $p_{\text {T}} ^{\gamma \gamma }$selection for masses below 400 GeV, and in Fig. 3b requiring $p_{\text {T}} ^{\gamma \gamma }$> 100 GeV for masses above 400 GeV, as well as for the search for non-resonant Higgs boson pair production. The fits are performed separately on the two distributions to search for resonant signals in both the low-mass and high-mass ranges. The observed data are found to be compatible with the sum of the expected SM backgrounds by performing a likelihood-ratio test [92]. The largest data excess has a local significance of 2.0 standard deviations at 400 GeV without the $p_{\text {T}} ^{\gamma \gamma }$selection. A modified frequentist method $\text {CL}_\text {s}$ [93] is used to calculate the 95% confidence-level (CL) exclusion limits with the asymptotic approximation [92]. Unfolding the SM Higgs boson branching fractions to $WW^*$ and $\gamma \gamma$ for the signal, the expected upper limit on the cross section for non-resonant Higgs boson pair production is 5.4 pb, while the observed limit is 7.7 pb, as shown in Table 6. The difference between the expected and observed limits is due to a slight excess of events in data. The expected upper limit on the cross section times the branching fraction of $X\rightarrow HH$ ranges from 17.6 to 4.4 pb, while the observed limit ranges from 40 to 6.1 pb, as a function of $m_X$ between 260 and 500 GeV, as shown in Fig. 4a.

Table 6 The 95% CL upper limits for the non-resonant production and the ratios of the limits to the SM cross-section value of $\sigma (pp\rightarrow HH)=33.4^{+2.4}_{-2.8}$ fb [17]. The $\pm 1\sigma$ and $\pm 2\sigma$ intervals around the median limit are also presented

Assuming the SM Higgs branching fractions of $B(H{\rightarrow }WW^*)=(21.52\pm 0.32)\%$ and $B(H{\rightarrow }{\gamma \gamma })=(0.227\pm 0.005)\%$ [17], the expected upper limit on the cross section for non-resonant production of $HH\rightarrow \gamma \gamma WW^*$ is 5.3 fb, while the observed limit is 7.5 fb, as shown in Table 6. The expected upper limit on the cross section for resonant production of $X\rightarrow HH\rightarrow \gamma \gamma WW^*$ ranges from 17.2 to 4.3 fb, while the observed limit ranges from 39.1 to 6.0 fb, as a function of $m_X$ between 260 and 500 GeV, as shown in Fig. 4b. The statistical uncertainty dominates in the final limits, while the impact of systematic uncertainties on these limits is only a few percent.

Fig. 4

$95\%$ CL expected (dashed line) and observed (solid line) limits on the resonant Higgs boson pair production cross section times the branching fraction of $X{\rightarrow }HH$ as a function of $m_X$ (a) with and b without assuming the SM branching fractions of $H{\rightarrow }WW^*$ and $H{\rightarrow }{\gamma \gamma }$. The $\text {CL}_\text {s}$ method and the asymptotic approximation are used. The $\pm 1\sigma$ and $\pm 2\sigma$ bands on the expected limit are also presented. To the right of the vertical dashed line at $m_{X} = 400$ GeV, the $p_{\text {T}} ^{\gamma \gamma } > 100$ GeV selection is applied in both plots, but not to the left

8 Conclusion

This paper presents searches for non-resonant and resonant Higgs boson pair production with a semileptonic $\gamma \gamma WW^*$ final state using 36.1 $\hbox {fb}^{-1}$ of $pp$ collision data collected at 13 TeV with the ATLAS detector at the LHC. No significant excess above the expected SM background is observed. A 95% confidence-level upper limit of 7.7 pb is set on the cross section for non-resonant production, while the expected limit is 5.4 pb, compared to the SM Higgs boson pair production cross section of 33.4 fb. The observed upper limit on the resonant production cross section times the branching fraction of $X\rightarrow HH$ ranges between 40 pb and 6.1 pb, while the expected limit ranges between 17.6 and 4.4 pb, for a hypothetical resonance with a mass in the range of 260–500 GeV. When deriving the limits above, the SM branching ratios of the $H{\rightarrow }{\gamma \gamma }$ and $H{\rightarrow }WW^*$ are assumed.