Search for Higgs boson pair production in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma \gamma WW^*$$\end{document}γγWW∗ channel using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$pp$$\end{document}pp collision data recorded at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}$$\end{document}s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$=13$$\end{document}=13 TeV with the ATLAS detector

Searches for non-resonant and resonant Higgs boson pair production are performed in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma \gamma WW^*$$\end{document}γγWW∗ channel with the final state of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma \gamma \ell \nu jj$$\end{document}γγℓνjj using 36.1 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {fb}^{-1}$$\end{document}fb-1 of proton–proton collision data recorded at a centre-of-mass energy of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=$$\end{document}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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma (pp\rightarrow X)\times B(X\rightarrow HH)$$\end{document}σ(pp→X)×B(X→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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H{\rightarrow }{\gamma \gamma }$$\end{document}H→γγ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H{\rightarrow }WW^*$$\end{document}H→WW∗ are assumed.


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 nonresonant and resonant production of Higgs boson pairs (H H) in the semileptonic γ γ W W * final state using 36.1 fb −1 of proton-proton ( pp) collision data recorded by the ATLAS detector at a centre-of-mass energy of √ s = 13 TeV. Previous searches for Higgs boson pair production have been performed by both the ATLAS and CMS experiments with data recorded at √ s = 8 TeV in the final states bbbb [8], bbγ γ [9,10], bbτ + τ − [11-13] and γ γ W W * [11], as well as e-mail: atlas.publications@cern.ch multi-lepton and multi-photon [14]. The pp collision data at √ s = 13 TeV have been analysed in order to search for Higgs boson pairs in the final states bbbb [15] and bbW W * [16]. No significant excess was observed compared to the SM prediction. However, it is important to explore the 13 TeV data in the channels that are not covered yet, such as the γ γ W W * 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 H H 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 tt H or H H H 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 heavyquark loop, where the cross section scales with the squared value of the tt H 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 (λ H H H ), 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 tt H , bbH or λ H H H 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 γ γ W W * final state (γ γ ν j j), i.e. with two photons, two jets, one charged lepton and a neutrino. This final state benefits from the large branching fraction of H →W W * [17], a characteristic signature from two photons and one lepton, as well as the excellent resolution of the diphoton invariant mass m γ γ , which provides good discrimination from a smooth continuum background composed of multi-photon and multi-jet SM processes. Given the expected sensitivity in 13 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.

The ATLAS detector
The ATLAS experiment [28] is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and nearly 4π coverage in solid angle. 1 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 |η| < 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 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 z axis along the beam pipe. The x axis points from the IP to the centre of the LHC ring, and the y axis points upwards. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). Angular distance is measured in units of 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 (|η| < 1.7). The endcap and forward regions are covered by LAr calorimeters for EM and hadronic energy measurements up to |η| = 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 twolevel trigger system is used to select events [30]. The firstlevel 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.

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 √ 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 ± 0.8 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 T > 35 GeV and the second with E 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.  [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 2 [36], and the parton distribution function (PDF) set CTEQ6L1 [37] was used.
• For non-resonant Higgs boson pair production, the inclusive cross sections 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-nextto-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.  ring 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.

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.
T is defined as the scalar sum of the transverse momenta ( p T ) of tracks with p T > 1 GeV within R = 0.2 of the photon candidate, excluding tracks from photon conversions. The selection on the track-based isolation variable is p iso Only photon candidates with |η| < 2.37 are considered, excluding the transition region between the barrel and endcap calorimeters (1.37 < |η| < 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 MS − q/ p 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 T or p 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 T of tracks lying within a cone of R = min (10 GeV/ p e(μ) T , 0.3) and the sum of E T of topological clusters of calorimeter cells within a cone of 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.057 p μ T + 95.57%) [72,73]. Finally, the electron candidates are required to have E T > 10 GeV and |η| < 2.47, excluding the transition region between the barrel and endcap calorimeters (1.37 < |η| < 1.52), whereas the muon candidates are required to have p T > 10 GeV and |η| < 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 |η| < 2.5 and p T > 25 GeV. In addition, a jet-vertex tagging algorithm (JVT) [81] is applied to jets with |η| < 2.4 and p 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-jets 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 btagging working point is selected to have an efficiency of 70% for a b-jet from tt 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 (lightflavour 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 R(e, γ ) < 0.4 are removed. Jets with R(jet, γ ) < 0.4 or R(jet, e) < 0.2 are removed. Electrons with R(e, jet) < 0.4 are removed. Muons with R(μ, γ ) < 0.4 or R(μ, jet) < 0.4 are removed.

Event selection
The events passing the diphoton trigger are required to contain at least two jets, no b-jet, and at least one charged lepton (e or μ, but including contributions from fully leptonic τ -lepton decays) in the final state. The two photon candidates with the leading (sub-leading) E T are required to satisfy E γ T /m γ γ > 0. 35 (0.25). The b-jet veto suppresses the tt H process. Furthermore, the transverse momentum of the diphoton system ( p γ γ T ) 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 γ γ T 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 GeV < m γ γ < 160 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 γ γ 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 γ γ distribution. The combination of two jets and at least one charged lepton would be consistent with H →W W * for the other Higgs boson. Its signature would in principle be enhanced by a missing transverse energy (E miss T ) requirement to indicate a neutrino, but a selection on E miss T was found not to produce any significant improvement in sensitivity, and so was not applied. The magnitude of E miss T [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.
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.

Signal and background estimation
A fit to the m γ γ 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 γ γ 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 χ 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  The theoretical uncertainties in the efficiency times acceptance ( × 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 lowp 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, α 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 →γ γ and H →W W * branching ratios (B) are considered as well. They are 2.1% and 1.5% [17], respectively.

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 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 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 ss | value reaches as large as 0.46 when not applying the p γ γ T selection, and 0.26 when applying the selection. The simulated background-only samples include the irreducible process of γ γ ν j j 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 γ γ T selection) and normalised according to the number of selected data events.

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 cal-ibration 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 γ events, and of electrons from Z → 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 γ γ 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 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 highp T ( p T > 300 GeV) jets and the other one is from extrapolating c-jets to τ -jets. The values for resonant production shown here assume m X = 260 GeV. Several theoretical uncertainties are reported for the cross section (σ ) and the combined efficiency and acceptance ( × A) 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.

Results
A fit to m γ γ is performed in the signal region to extract the signal yield. The statistical model is constructed with a likelihood function: • i stands for the event index,  G(0|θ, 1) is the pdf of a Gaussian distribution used to constrain the nuisance parameters θ 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. 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 γ γ T selection for masses below 400 GeV, and in Fig. 3b requiring p γ γ T > 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 γ γ T selection. A modified frequentist method CL s [93] is used to calculate the 95% confidencelevel (CL) exclusion limits with the asymptotic approximation [92]. Unfolding the SM Higgs boson branching fractions to W W * and γ γ 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 → H H 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.
Assuming the SM Higgs branching fractions of B (H →W W * ) = (21.52 ± 0.32)% and B(H →γ γ ) = (0.227 ± 0.005)% [17], the expected upper limit on the cross section for non-resonant production of H H → γ γ W W * 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 → H H → γ γ W W * 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  The CL s method and the asymptotic approximation are used. The ±1σ and ±2σ bands on the expected limit are also presented. To the right of the vertical dashed line at m X = 400 GeV, the p γ γ T > 100 GeV selection is applied in both plots, but not to the left limits, while the impact of systematic uncertainties on these limits is only a few percent.

Conclusion
This paper presents searches for non-resonant and resonant Higgs boson pair production with a semileptonic γ γ W W * final state using 36.1 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 → H H 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 →γ γ and H →W W * are assumed.