Search for Higgs boson pair production in the $b\bar{b}WW^{*}$ decay mode at $\sqrt{s}=13$ TeV with the ATLAS detector

A search for Higgs boson pair production in the $b\bar{b}WW^{*}$ decay mode is performed in the $b\bar{b}l\nu qq$ final state using 36.1 fb$^{-1}$ of proton-proton collision data at a centre-of-mass energy of 13 TeV recorded with the ATLAS detector at the Large Hadron Collider. No evidence of events beyond the background expectation is found. Upper limits on the non-resonant $pp\to HH$ production cross section of 10 pb and on the resonant production cross section as a function of the $HH$ invariant mass are obtained. Resonant production limits are set for scalar and spin-2 graviton hypotheses in the mass range 500 to 3000 GeV.


Introduction
The Higgs boson (H) is an essential part of the Standard Model (SM) and it has a crucial role in the electroweak symmetry breaking (EWSB) mechanism [1-6]. In this mechanism, an SU(2) doublet bosonic scalar field is subject to a potential energy term whose shape allows the doublet field to acquire a vacuum expectation value that breaks the SU(2) symmetry and produces the Higgs boson and its potential energy term. This potential is the last piece of the SM Lagrangian which is yet to be directly tested.
The shape of the Higgs boson potential in the SM can be expressed as a function of the Fermi coupling constant G F and the Higgs boson mass m H . A direct phenomenological prediction of the SM due to the potential is the interaction of the Higgs boson with itself at tree level (self-interaction), which can be probed by studying di-Higgs boson production in proton-proton collisions, as illustrated in Figure 1(a). The self-interaction diagram together with the quark-loop contributions, primarily via the top-Higgs Yukawa coupling, Figure 1(b), are the leading-order Feynman diagrams for Higgs boson pair production. The SM cross section for pp → HH is extremely small, e.g. 33.4 fb at 13 TeV [7]. Physics beyond the SM can manifest in the increased production with respect to the SM predictions of the non-resonant HH final state or in the resonant production of particles that decay into a pair of SM Higgs bosons. The analysis presented here is potentially sensitive to cases where the decaying particle is a scalar, as in the MSSM [8] and 2HDM models [9], or a spin-2 graviton, as in Randall-Sundrum models [10]. The signals under study are non-resonant HH production with event kinematics predicted by the SM and resonant HH production with event kinematics consistent with the decays of heavy spin-0 or spin-2 resonances.
Results at √ s = 13 TeV were published by the ATLAS Collaboration in the 4b [17], bbτ + τ − [18], bbγγ[19] and WW γγ[20] decay mode and by CMS in the 4b [21] , bbτ + τ − [22], bbγγ [23] and in the bbWW * channel using the dileptonic WW * decay mode [24]. Given the low expected yield for SM HH nonresonant production, it is of great importance to understand the sensitivity for the observation of the Higgs boson pair production in all possible decay channels, including bbWW * , which will improve projections for future high-luminosity and high-energy colliders. This paper reports results of a search for Higgs boson pair production where one Higgs boson decays via H → bb, and the other decays via H → WW * . The H → WW * branching fraction is the second largest after H → bb, so the bbWW * final state can be sensitive to HH production if the signal can be well separated from the dominant tt background. The WW * system decays into νqq (where is either an electron or a muon), and the small contamination from leptonic τ decays is not explicitly vetoed in the analysis. Figure 2 shows a schematic diagram of resonant production of the Higgs boson pair with the subsequent decays H → WW * and H → bb.
Two complementary techniques are used to reconstruct the Higgs boson candidate that decays into two b-quarks. Both techniques use the anti-k t jet algorithm [25] but with different radius parameters. The first technique employs jets with radius parameter R = 0.4 and it is used when each b-quark from the H → bb decay can be reconstructed as a distinct b-jet. The second technique is used when this is not possible, due to the large boost of the b-quark pair. In this case the Higgs boson candidate is identified as a single anti-k t jet with radius parameter R = 1.0. The analysis using the first technique is referred to as the "resolved" analysis and that using the second technique is referred to as the "boosted" analysis. In both analyses, the jets from the hadronically decaying W boson are reconstructed as anti-k t jets with radius parameter R = 0.4. The resonant HH search is performed using either the resolved or boosted analysis method depending on which is the most sensitve to the particular model and HH mass being tested, in contrast to the non-resonant search which uses only the resolved analysis method. The dominant background in the bbWW * final state is tt production, with smaller contributions from W bosons produced in association with jets (W+jets) and multijet events in which a jet is misidentified as a lepton. The analysis defines one signal region for each signal hypothesis and, in order to avoid biases in the analysis selection, the analysis procedures and the event selection were optimised without reference to data in the signal regions.

Data and simulation samples
The ATLAS detector [26] is a general-purpose particle detector at the Large Hadron Collider optimised to discover and measure a broad range of physics processes. It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer incorporating three large superconducting toroid magnets.1 The dataset used in this analysis corresponds to an integrated luminosity of 36.1 fb −1 (3.2 fb −1 from 2015 and 32.9 fb −1 from 2016) recorded by single-electron or single-muon triggers. The single-lepton trigger efficiency ranges from 75% to 90% (75% to 80%) for electrons (muons) depending on the signal mass, for selected lepton candidates above p T thresholds defined in Section 4.1. Samples of simulated signal and background events were used to design the event selection and estimate the signal acceptance and the background yields from various SM processes.
When searching for a new resonance (denoted by X in the following), specific simulation models must be employed. Therefore, the spin-0 states were treated as narrow heavy neutral Higgs bosons, while the spin-2 states were modelled as narrow Randall and the UE-EE-5-CTEQ6L1 tune. The simulation produced the Higgs boson pair through gluon-gluon fusion using an effective field theory approach to take into account the finite value of the top-quark mass m t [35]. Events were first generated with an effective Lagrangian in the infinite top-quark mass approximation, and then reweighted with form factors that take into account the finite mass of the top quark.
The non-resonant signal samples were simulated with M 5_ MC@NLO + H ++ using the CT10 PDF set; and the same approach for the inclusion of finite m t effects was used [36]. In addition, scale factors dependent on the HH invariant mass m H H at generator level were applied to match the MC m H H distribution with an NLO calculation that computes exact finite m t contributions [37]. All signal samples were generated with 100% of Higgs boson pairs decaying into bbWW * , and the samples were then normalised assuming B(H → WW * ) = 0. 22  The tt background samples were generated with P -B v2 [41] using the CT10 PDF set. P -B v2 was interfaced to P 6.428 [42] for parton showers, using the P 2012 [43] tune with the CTEQ6L1 [44] set of PDFs for the underlying-event description. E G v1.2.0 [45] was used to simulate the bottom and charm hadron decays. The mass of the top quark was set to m t = 172.5 GeV. At least one top quark in the tt event was required to decay into a final state with a lepton. For the 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). The angular distance is measured in units of ∆R ≡ (∆η) 2 + (∆φ) 2 .
tt sample the parameter H , used to regulate the high-p T gluon emission in P , was set to m t , giving good modelling of the high-p T region [46]. The interference between the tt background and the signal is extremely small due to the small width of the Higgs boson (Γ H ∼4 MeV) and it has been neglected in this analysis. The tt cross section is calculated to next-to-next-to-leading order in QCD including resummation of soft gluon contributions at next-to-next-to-leading-logarithm (NNLL) accuracy using T ++ 2.0 [47].
Single-top-quark events in the Wt-, s-and t-channels were generated using P -B v1 [48,49]. The overall normalisation of single-top-quark production in each channel was rescaled according to its approximate NNLO cross section [50][51][52].
The effect of multiple pp interactions in the same and neighbouring bunch crossings (pile-up) was included by overlaying minimum-bias collisions, simulated with P 8.186, on each generated signal and background event. The interval between proton bunches was 25 ns in all of the data analysed. The number of overlaid collisions was such that the distribution of the number of interactions per pp bunch crossing in the simulation matches that observed in the data: on average 14 interactions per bunch crossing in 2015 and 23.5 interactions per bunch crossing in 2016. The generated samples were processed through a G 4-based detector simulation [53,54] with the standard ATLAS reconstruction software used for collision data.

Object reconstruction
In the present work an "object" is defined to be a reconstructed jet, electron, or muon. Electrons are required to pass the "TightLH" selection as described in Ref. [55,56], have p T > 27 GeVand be within |η| < 2.47, excluding the transition region between the barrel and endcaps in the LAr calorimeter (1.37 < |η| < 1.52). In addition, the electron is required to be isolated. In order to calculate the isolation variable, the p T of the tracks in a cone of ∆R around the lepton track is summed ( p T ), where ∆R = min(10 GeV/p e T , 0.2) and p e T is the electron transverse momentum. The ratio p T /p e T (isolation variable) is required to be less than 0.06.
Muons are reconstructed as described in Ref. [57] and required to pass the "Medium" identification criterion and have |η| < 2.5. The muon isolation variables are similar to the electron isolation variables with the only difference being that the maximum cone size is ∆R = 0.3 rather than 0.2.
Jets are reconstructed using the anti-k t algorithm [25] with a radius parameter of 0.4, and are required to have p T > 20 GeV and |η| < 2.5. Suppression of jets likely to have originated from pile-up interactions is achieved using a boosted decision tree in an algorithm that has an efficiency of 90% for jets with p T < 50 GeV and |η| < 2.5 [58].
Large-R jets are reconstructed using the anti-k t algorithm with a radius parameter of 1.0 and are trimmed to reduce pile-up contributions to the jet, as described in Ref. [59]. The jet mass (m J ) resolution is improved at high momentum by using tracking in addition to calorimeter information [60]. Large-R jets are required to have p T > 250 GeV, m J > 30 GeV and |η| < 2.0. The identification of large-R jets consistent with boosted Higgs boson decays uses jets built from ID tracks (referred to as track-jets) to identify the b-jets within the large-R jets. The track-jets are built with the anti-k t algorithm with R = 0.2 [61]. They are required to have p T > 10 GeV, |η| < 2.5, and are matched to the large-R jets with a ghost-association algorithm [62].
The jet-flavour tagging algorithm [63] is used to select signal events and to suppress multijet, W+jets, Z+jets and diboson backgrounds. The jets containing b-hadrons are called b-jets in this work. The jet-flavour tagging algorithm parameters were chosen such that the b-tagging efficiency is 85% for jets with p T of at least 20 GeV as determined in simulated inclusive tt events [63] . At this efficiency, for jets with a p T distribution similar to that originating from jets in tt events, the charm-quark component is suppressed by a factor of 3.1 while the light-quark component is suppressed by a factor of 34. Jets that are not tagged as b-jets are collectively referred to as "light jets".
The calorimeter-based missing transverse momentum with magnitude E miss T is calculated as the negative vectorial sum of the transverse momenta of all calibrated selected objects, such as electrons and jets, and is corrected to take into account the transverse momentum of muons. Tracks with p track T > 500 MeV, compatible with the primary vertex but not matched to any reconstructed object, are included in the reconstruction to take into account the soft-radiation component that does not get clustered into any hard object [64].
To avoid double-counting, overlapping objects are removed from the analysis according to the following procedure. Muons sharing their track with an electron are removed if they are calorimeter-tagged. Otherwise, the electron is removed. Jets overlapping with electrons within an angular distance ∆R = 0.2 are removed. Jets overlapping with muons within ∆R = 0.2 and having less than three tracks or carrying less than 50% of the muon p T are removed. Electrons overlapping with remaining jets within ∆R = min(0.4, 0.04 + 10 GeV/p e T ) are removed. Muons overlapping with remaining jets within ∆R = min(0.4, 0.04 + 10 GeV/p µ T ) are removed. s 4 Resolved analysis

Resolved analysis: event selection
At lowest order in QCD the final-state particles consist of one charged lepton, one neutrino, and jets of colourless hadrons from four quarks, two being b-quarks. Therefore, the corresponding detector signature is one charged lepton (e/µ), large E miss T , and four or more jets. Two of these jets are b-tagged jets from the Higgs boson decay, and two jets are non-b-tagged jets from the hadronic W boson decay.
The data used in the analysis were recorded by several single-electron or single-muon triggers in 2015 and 2016. In 2015, the electron (muon) trigger required a p T > 24 (20) GeV electron (muon) candidate. Because of a higher instantaneous luminosity, in 2016 the electron trigger required a p T > 26 GeV electron candidate, while muons were triggered using a p T threshold of 24 GeV at the beginning of data taking, and 26 GeV for the rest of the year. In both 2015 and 2016, a threshold of p T > 27 GeV was applied offline on the selected lepton candidate.
The analysis selects events that contain at least one reconstructed electron or muon matching a trigger lepton candidate. In order to ensure that the leptons originate from the interaction point, requirements on the transverse (d 0 ) and longitudinal (z 0 ) impact parameters of the leptons relative to the primary vertex are imposed. In particular, defining σ d 0 as the uncertainty in the measured d 0 and θ as the angle of the track relative to the beam axis, the requirements |d 0 |/σ d 0 < 2 and |z 0 sin θ| < 0.5 mm are applied. The requirement on |d 0 |/σ d 0 is relaxed to define control regions in order to estimate the multijet background. The highest p T lepton is then retained as the analysis lepton. Events are required to have exactly two b-tagged jets, which form the Higgs boson candidate. Since events are accepted if they contain two or more light jets, in events with more than two light jets, the three leading jets are considered, and the pair with the lowest ∆R between them is selected as the W boson candidate. From MC simulation it was found that, when the light quarks from the W boson are matched to reconstructed jets by requiring that the ∆R between the jet and the quark is less than 0.3, this procedure yields the correct jet assignment in 70% of the cases.
The event kinematics of the H → WW * → νqq topology can be fully reconstructed. Among all fourmomenta of the final-state particles, only the component of the neutrino momentum along the beam axis, referred to as longitudinal momentum (p z ) in the following, is unknown while its transverse momentum is assumed to be the E miss T . The neutrino longitudinal momentum is computed by solving a quadratic equation in p z , employing the four-momenta of the lepton and the hadronic W boson, the E miss T , and the m H = 125 GeV constraint on the WW * system. No W boson mass constraint is applied to either the hadronic or the leptonic W boson decay, allowing either W boson to be off-shell. Whenever two real solutions are obtained, the ν candidate with the smallest ∆R relative to the lepton direction is retained. Studies performed by matching the ν candidate with the MC generator-level neutrino show that this procedure finds the correct solution for the neutrino p z in 60% (75%) of cases for a resonant signal of mass 700 (3000) GeV. If two complex solutions are found, only the real part of the solutions is retained. With the neutrino longitudinal momentum computed, the di-Higgs invariant mass can be fully reconstructed and employed to discriminate against backgrounds.
Kinematic selections are used to suppress the tt background relative to the signal. The tt events are typically characterised by two b-jets and two W bosons such that the ∆R separation between the b-jets is large, and similarly the ∆R separation between the W bosons is also large. In contrast, in particular when the invariant mass of the heavy resonance is large, the signal is characterised by two b-jets and two W bosons which are closer in ∆R in signal events with respect to the tt background events. Moreover, for the signal the two b-jets have an invariant mass equal to m H , while this is not the case for the tt background, where a much broader distribution is expected. The symbols of the kinematic variables that discriminate between signal and background are listed in Table 1.
The selection requirements on the kinematic variables defining the signal region were chosen to maximise the expected sensitivity to various signals. The optimisation was performed for a spin-0 signal considering resonance masses (m X ) from 500 GeV to 3000 GeV in steps of 100 GeV. The same selection was used for the spin-2 signal models while SM Higgs pair production was used to optimise the non-resonant analysis. Below 500 GeV the top-quark background increases significantly, and hence rapidly reduces sensitivity.
105-135 105-135 105-135 105-135 The selection criteria define four sets of requirements, referred as non-res, m500, low-mass and high-mass in the following. They are shown in Table 2. The non-res and m500 selections are exclusively used for non-resonant signal and resonant signal with mass 500 GeV respectively. The low-mass selection is used for signal masses from 600 to 1300 GeV, while the high-mass selection is used for signals with masses between 1400 and 3000 GeV. In addition, requirements are placed on the reconstructed di-Higgs invariant mass m H H as a function of the signal resonance mass m X , as shown in Table 3. The resolution of the reconstructed m H H ranges from 6% at 500 GeV to 10% at 3000 GeV.

Resolved analysis: background determination
In this analysis the presence of a signal is indicated by an excess of events over the SM prediction for the background yield in the signal regions, so it is of great importance to properly estimate the amount of background in those regions. The dominant background is the tt process. Dedicated control regions are used to normalise and validate the estimate of this background. The tt normalisation is performed using three data control regions, one for the non-res, a second one for the m500 and low-mass, and a third one for the high-mass selection. These control regions are obtained by selecting events outside the m bb window [100, 140] GeV and applying only the E miss T , m WW * (where applicable) and p bb T requirements shown in Table 2 for the respective selections.
In all regions, the event yields of W/Z+jets, single-top-quark and diboson events are modelled using simulated events and normalised to the expected SM cross sections. The multijet component of the background originates from events where either a jet is incorrectly identified as a lepton, or a non-prompt lepton is produced in heavy-flavour decays, or from photon conversions. It is characterised by low E miss T and high |d 0 |/σ d 0 values of the lepton. The multijet background makes a significant contamination in the top control regions. Therefore, this background is estimated in each top control region and signal region using a data-driven two-dimensional sideband method, labelled the ABCD method, that uses three additional regions denoted in the following by B, C and D. The region of interest, signal or control region, is indicated by A.
The B, C and D regions are defined in the following way: while N A ,N B ,N C and N D indicate the number of events in the A,B,C and D regions respectively. In the absence of correlations between the E miss T and |d 0 |/σ d 0 variables, the relation N A = N C N B /N D holds, while in practice a correlation among variables results in a correction factor F to be applied to the computed ratio N corrected The correction factor F is estimated from data at an early stage of the analysis selection once a veto on the signal candidates is applied by inverting the requirement on the m bb variable. It is computed using the relation F = N A N D /(N C N B ). Systematic uncertainties in F are described in Section 4.3. In order to reduce statistical uncertainties in the computation, the shape of the m bb distribution is derived at an earlier stage of the selection sequence, after applying the m WW * < 130 GeV and p bb T > 210 GeV requirements for the non-res, m500 and low-mass analyses and the p bb T > 350 GeV and p WW * T > 250 GeV requirements for the high-mass analysis. It was verified that subsequent requirements do not affect the m bb shape, which can therefore be used at the end of the selection sequence. Table 4 summarises the numbers of observed and estimated events in the three top-quark control regions. The event yields in the control regions are used as input to the statistical analysis. Major contamination in the tt control regions comes from multijet and W+jets backgrounds; as a result the tt purity ranges from 52% to 58%.  The modelling of the background was checked at all selection stages and, in general, shows good agreement with data. Figure 3 shows the m T distribution of the leptonic W boson candidate in the three top control regions. The m T variable is defined as: where ∆φ is the azimuthal angle between p T and E miss T . The multijet background populates the low values of the m T distribution, so any mis-modelling of the multijet background would be clearly visible in the m T distribution.   The m bb distribution in the resolved analysis for the non-res and m500 selections at the end of the selection sequence, before applying the m bb requirement. The signals shown are from SM non-resonant HH production scaled up by a factor of 300 (left) and from a scalar resonance with mass 500 GeV scaled to the expected upper-limit cross section reported in Section 6 (right). The lower panel shows the fractional difference between data and the total expected background with the corresponding statistical and total uncertainty.

Resolved analysis: systematic uncertainties
The main systematic uncertainties in the background estimate arise from the potential mis-modelling of background components. For tt background, MC simulation is used to derive the acceptances in all analysis regions, while the normalisation is taken from the top control region and applied in the signal regions. Therefore, the acceptance ratio between signal and control regions is affected by theoretical uncertainties in the simulated tt sample. These uncertainties are estimated by considering five sources: the matrix element generator used for the tt simulation and the matching scheme used to match the NLO matrix element with the parton shower, the parton shower modelling, the initial-state (Initial State Radiation, ISR) and final-state (Final State Radiation, FSR) gluon emission modelling, the dependence on the choice of the PDF set and the dependence on the renormalisation and factorisation scales. Matrix element generator and matching systematic uncertainties are computed by comparing samples generated by aMC@NLO [29] and P , both interfaced with H ++ for showering and fragmentation. Parton shower systematic uncertainties are computed by comparing samples generated using P +P 6 and P +H ++. Initial-state and final-state radiation systematic uncertainties are computed by varying the generator parameters from their nominal values to increase or decrease the amount of radiation. The PDF uncertainties are computed using the eignevenctors of the CT10 PDF set. Uncertainties due to missing higher-order corrections, labelled scale uncertainties, are computed by independently scaling the renormalisation and factorisation scales in aMC@NLO+H ++ by a factor of two, while keeping the renormalisation/factorisation scaling ratio between 1/2 and 2. These systematic uncertainties are summarised in Table 5.
Uncertainties in the modelling of W+jets background are computed in each signal region (SR) and top control region (CR). Three sources of uncertainty are considered: scale variation, PDF set variation and generator modelling uncertainties. Scale uncertainties are computed by scaling the nominal renormalisation and factorisation scales by a factor of two. PDF uncertainties are computed using the NNPDF [39] error set, while generator modelling uncertainties are obtained by comparing the nominal S generated sample with a sample generated with A [65] and showered with P 6 [42]. The values obtained in each region are summarised in Table 6.
For the data-driven multijet background, three sources of uncertainty are identified. The non-closure correction term F is computed using data at an early stage of the selection sequence, where contamination by the signal can be considered negligible. Its difference from the value obtained using a simulated multijet event sample is 40% and is assigned as an uncertainty in the multijet estimation. The F value can be affected by the analysis selection requirements. A systematic uncertainty (extrapolation uncertainty) is added by comparing the maximum variation among the F values evaluated after each selection requirement. Finally,  Shower 40  40  40  40  20  20  PDF  30  7  40  10  30  20  Scale  20  30  20  30  30  30 the uncertainty due to the dependence of the F value on lepton flavour (flavour uncertainty) is computed as the maximum difference between the nominal F value and the F value calculated for electrons and muons separately. The extrapolation (flavour) uncertainty is found to be 16% (9%) for the non-res selection, 32% (9%) for the m500 and low-mass resonant selections, and 45% (6%) for the high-mass resonant selection.
Single-top-quark production is one of the smaller backgrounds in this analysis. Theoretical cross-section uncertainties vary from 5% for associated Wt production to 4% for s-and t-channel single-top production.
The largest of these is conservatively assigned to all single-top production modes. Further modelling systematic uncertainties are calculated by employing the difference between the nominal sample using the Diagram Removal scheme described in Ref.
[66] and a sample using the Diagram Subtraction scheme for the dominant single-top production mode, Wt. The uncertainties are 50%, for the non-res, m500 and low-mass analyses, and 80% for the high-mass analysis.
Systematic uncertainties in the signal acceptance are computed by varying the renormalisation and factorisation scales with a variation of up to a factor of two, and using the same procedure as for the tt background. PDF uncertainties are computed using PDF4LHC15_30 [67] PDF sets, which include the envelope of three PDF sets, namely CT14, MMHT14, NNPDF3.0. The resulting uncertainties are less than 1.1% for the scale and less than 1.3% for the PDFs. Parton shower uncertainties are computed by comparing the H ++ showering with that of P 8, and this results in less than 2% uncertainty.

Boosted analysis: event selection
As in the resolved analysis, data used in the boosted analysis were recorded by single-lepton triggers, and only events that contain at least one reconstructed electron or muon matching the trigger lepton candidate are analysed. Requirements on p T , |d 0 |/σ d 0 and |z 0 sin θ| of the lepton tracks are also the same as in the resolved analysis.
Events are required to have at least one large-R jet with an angular distance ∆R > 1.0 from the reconstructed lepton. The highest-p T large-R jet is identified as the H → bb candidate. The large-R jet mass is required to be between 30 GeV and 300 GeV. In order to reconstruct the H → WW * system, events with at least two small-R jets with an angular distance ∆R > 1.4 from the H → bb candidate are selected. The hadronically and leptonically decaying W bosons are then reconstructed following the same algorithm as in the resolved analysis. In order to reduce the tt background, events are rejected if they contain any small-R jet passing the b-tagging requirement.
Signal regions (SR) are defined with at least two associated track-jets within the large-R jet and requiring that the two highest-p T track-jets are also b-tagged. The large-R jet mass must be between 90 GeV and 140 GeV. An additional requirement of E miss T > 50 GeV is imposed to reject multijet backgrounds.
In order to assess the modelling of the dominant tt background, a validation region (VR) is defined outside the large-R jet signal region mass window and labelled top VR. Any event with a large-R jet mass m Large-R jet < 90 GeV or m Large-R jet > 140 GeV falls in the top VR. By construction, the top VR is orthogonal to the SR.

Boosted analysis: background determination
In the boosted analysis the presence of a signal is indicated by an excess of events above the SM prediction of the background m H H distribution at the end of the event selection. Similarly to the resolved analysis, the tt process is the dominant background. Therefore, a dedicated validation region is used to check its modelling as defined in Section 5.1. The event yields from tt, W/Z+jets, single-top-quark and diboson processes in the signal region and the top VR are modelled using simulation and normalised to the expected SM cross section described in Section 2.
The multijet component of the background is estimated using the data-driven method as in the resolved analysis. In the boosted analysis a higher requirement on E miss T (E miss T > 50 GeV) is applied, while the cut on |d 0 |/σ d 0 is the same. For the boosted analysis, the correlation between |d 0 |/σ d 0 and E miss T is estimated in multiple MC background samples and also in data, and it is found to be negligible. Hence, the multijet yield in region A can be estimated using the relation N A = N C N B /N D . The multijet estimation is performed separately for the muon and the electron channel. The N B /N D ratio is calculated inclusively in the large-R jet mass distribution. The m H H distribution of the multijet background is estimated by subtracting the prompt-lepton MC backgrounds from the data in the 1-tag region, where the 1-tag region is defined as the region where all selections are applied except that the large-R jet is required to have only one track-jet tagged as a b-jet.
The modelling of the background is checked in the top VR. Table 7 reports the numbers of observed and predicted background events in the top VR, showing good agreement between the two. In order to check   The signal distribution is negligible in the left plot, while in the right plot it has been scaled to the expected upper-limit cross section reported in Section 6. The lower panel shows the fractional difference between data and the total expected background with the corresponding statistical and total uncertainty. the validity of the multijet background determination, the m T distribution is shown in Figure 6. This variable is particularly sensitive to the multijet background contamination. Additionally, the m Large-R jet variable used to define the signal region and the top VR is shown in the same figure. The data and predicted background agree well, which builds confidence in the estimated efficiency of the m Large-R jet requirement for signal and background.

Boosted analysis: systematic uncertainties
The evaluation of detector modelling uncertainties in the boosted analysis follows the same approach as in the resolved analysis. The significant additions to those described in Section 4.3 are the uncertainties related to the large-R jets. The large-R jet energy resolution and scale, and jet mass resolution and scale uncertainties are derived in situ from 8 TeV pp collision data, taking into account MC simulation extrapolations for the different detector and beam conditions present in 8 and 13 TeV data-taking periods [71].
The uncertainty in the b-tagging efficiency for track-jets is evaluated with the same method used for resolved calorimeter jets. The impact of these uncertainties on the final fit are shown in Table 13.
All SM backgrounds, except multijet, are modelled using MC simulation. Therefore, predicted yields in both the signal and the top validation regions are affected by theoretical uncertainties. These uncertainties are computed following the same procedure as in the resolved analysis for tt, W/Z+jets, single-top-quark and diboson backgrounds. For the tt background in the signal region, the uncertainties are summarised in Table 8. The uncertainties on single top quark production range from 20% for ISR/FSR to 70% stemming from the difference between the diagram removal and diagram subtraction schemes. Uncertainties in the modelling of W/Z+jets background range from 10% stemming from PDF uncertainties to 45% stemming from scale uncertainties. Diboson processes have a negligible impact on the total background. For the normalisation of the multijet background predicted in region A (See Section 5.2), several sources of uncertainty are considered. The uncertainties in the normalisation of tt and W/Z+jets in regions B, C and D contribute a systematic uncertainty of 25% and 30% respectively. The relative difference between the large-R jet mass acceptance in the 1-tag region C and in the 2-tag region C accounts for 15%. The propagation of the statistical uncertainty in the multijet yield in region C and the uncertainty in the N B /N D ratio contribute about 23%. The propagation of detector modelling systematic uncertainties, including the modelling uncertainty of the |d 0 |/σ d 0 requirement and of the MC backgrounds with prompt leptons subtracted from data in regions B, D and C, contribute about 45%. As an additional check on the prediction of the multijet yield with the ABCD method, a conditional background-only likelihood fit of the large-R jet mass distribution is performed in the VR. The difference between the multijet yield estimated with this method and the ABCD prediction is assigned as an uncertainty. This error accounts for 23% of the total uncertainty in the multijet estimation. All different sources of uncertainty are treated as independent and added in quadrature for the final uncertainty of 80% in the multijet normalisation.
For the simulated backgrounds, the systematic uncertainty in the m H H distribution shape is determined by comparing the nominal MC sample with the corresponding alternative (variation) MC samples described in Section 4.3. The shape systematic uncertainty is determined by fitting a first-order polynomial to the ratio of the variation m H H distribution to the nominal m H H distribution, while keeping the same normalisation. For the data-driven multijet background, the uncertainty in the m H H distribution shape is determined by comparing the shapes in the 2-tag and 1-tag C regions.
Systematic uncertainties in the signal acceptance are computed following the same algorithm as the resolved analysis. The resulting uncertainties are less than 0.5% for uncertainties due to missing higher-  order corrections (labelled scale), less than 0.5% for those due to PDFs, and approximately 2% (5%) in the lower (higher) mass range for those due to the parton shower.

Resolved analysis
After applying the requirements listed in Table 2, the invariant mass of the HH system (m H H ) is distributed as shown in Figures 7 and 8. Data are generally in good agreement with the expected background predictions within the total uncertainty. The signal m H H distribution is shown in the figure for the nonresonant, the scalar resonance, and the two graviton hypotheses with c = 1.0 and c = 2.0. The scalar samples are simulated in the narrow-width approximation, so the reconstructed width is exclusively due to the detector resolution. The same holds for graviton samples with c = 1.0, while c = 2.0 graviton samples have a significant intrinsic width that leads to a loss of sensitivity.
The m H H distribution is sampled with resonance-mass-dependent m H H requirements as reported in Table  3. The numbers of events in the signal and control regions (the tt control region and the C region of the multijet estimation procedure) are simultaneously fit using a maximum-likelihood approach. The fit includes six contributions: signal, W+jets, Z+jets, tt, single-top-quark production, diboson and multijet. The tt and multijet normalisations are free to float, the C region of the ABCD method being directly used in the fit, while the diboson, W+jets and Z+jets backgrounds are constrained to the expected SM cross sections within their uncertainties.
The fit is performed after combining the electron and muon channels. Statistical uncertainties due to the limited sample sizes of the simulated background processes are taken into account in the fit by means of nuisance parameters, which are parameterised by Poisson priors. Systematic uncertainties are taken into  account as nuisance parameters with Gaussian constraints. For each source of systematic uncertainty, the correlations across bins and between different kinematic regions, as well as those between signal and background, are taken into account. Table 9 shows the post-fit number of predicted backgrounds, observed data, and the signal events normalised to the expected upper limit cross sections.
No significant excess over the expectation is observed and the results are used to evaluate an upper limit at the 95% confidence level (CL) on the production cross section times the branching fraction for the signal hypotheses under consideration. The exclusion limits are calculated with a modified frequentist method [72], also known as CL s , and the profile-likelihood test statistic [73]. None of the considered systematic uncertainties is significantly constrained or pulled in the likelihood fit.
The branching fraction B(HH → bbWW * ) = 2 × B(H → bb) × B(H → WW * ) = 0.248 is used to obtain the following observed (expected) limit on the HH production cross section at 95% CL: σ(pp → HH) < 10 10 +4 −3 pb, which corresponds to 300 (300 +100 −80 ) times the SM predicted cross section. Including only the statistical uncertainty, the expected upper limit for the non-resonant production is 190 times the SM prediction. This result, when compared with other HH decay channels, is not competitive. This is mainly due to the similarity of the reconstructed m H H spectrum between the non-resonant SM signal and the tt background that makes the separation between the two processes difficult. Table 9: Data event yields, and post-fit signal and background event yields in the final signal region for the nonresonant analysis and the resonant analysis in the 500-3000 GeV mass range. The errors shown are the MC statistical and systematic uncertainties described in Section 4.3. The yields are shown for three signal models: a scalar (S) and two Randall-Sundrum gravitons with c = 1.0 and c = 2.0 (G * KK ). Signal event yields are normalised to the expected upper-limit cross section. 18.4 ± 1.5 19.7 ± 1.6 18.2 ± 1.5 20 ± 8 28 900
The analysis is most sensitive for a mass value of 1300 GeV with an expected upper limit of 0.35 pb on σ(pp → HH). At this mass the observed exclusion limit is 0.2 pb. In both the non-resonant and resonant cases, the impact of the systematic uncertainties is observed to be large. In order to quantify the impact of the systematic uncertainties, a fit is performed where the estimated signal yield, normalised to an arbitrary cross-section value, is multiplied by a scaling factor α sig , which is treated as the parameter of interest in the fit. The fit is performed using pseudo-data and the contribution to the uncertainty in α sig from several sources is determined. The contribution of the statistical uncertainty to the total uncertainty in α sig , shown in Table 10, is decomposed into signal region statistics, top CR statistics and multijet CR statistics. The contribution of the systematic uncertainties to the total uncertainty is decomposed into the dominant components and shown in Table 11. The dominant systematic uncertainties vary across the mass range, but some of the most relevant ones are due to tt modelling, b-tagging systematic uncertainties, and those related to jet measurements.   Table 11: Systematic contributions (in percentage) to the total error in the scaling factor α sig for the non-resonant signal and three scalar-signal mass hypotheses, 500 GeV, 1000 GeV and 2000 GeV, in the resolved analysis. The first column quotes the source of the systematic uncertainty. The " − " symbol indicates that the specified source is negligible. The contribution is obtained by calculating the difference in quadrature between the total error in α sig and that obtained by setting constant the nuisance parameter(s) relative to the contribution(s) under study.

Boosted analysis
The boosted analysis applies the selection criteria described in Section 5.1. After applying the large-R jet mass requirement 90 < m Large-R jet < 140 GeV, the m H H distribution is reconstructed and its shape is fit to data using MC signal and background templates. The distribution is fit using 17 bins, with almost uniform width except at low and high m H H , where the bin width is modified in order to have a MC statistical uncertainty smaller than 20%. All backgrounds, except multijet, are simulated using MC generators and normalised using the cross section of the simulated process. The multijet background is estimated using the ABCD method, and its normalisation obtained from this method is kept fixed in the fit. The bias due to possible signal contamination in the ABCD regions was studied and found to have negligible effect on the result. The integral of the m H H distribution for the boosted analysis is shown in Table 12.
Systematic uncertainties affecting the m H H shape are parameterised as linear functions of m H H , and the function parameters are treated as nuisance parameters in the fit. Statistical uncertainties due to the limited sample sizes of the simulated background processes are taken into account in the fit by means of further nuisance parameters, which are parameterised by Poisson priors.
The systematic uncertainties included in the fit are described in Section 5.3. The contribution of the systematic uncertainties to the total uncertainty is decomposed into the dominant components and summarised in Table 13. The most relevant systematic uncertainties are due to the limited size of the MC samples, the tt modelling and the b-tagging systematic uncertainties.  Figure 11 shows the observed and the expected upper limit on the production cross section of the scalar S and the two graviton (c = 1.0 and c = 2.0) G * KK particles.

Summary
The boosted and resolved analyses have non-trivial event overlap, so their limits are jointly represented through a priority scheme where the two analyses are presented in the mass range where they are most sensitive. The boosted analysis was studied for masses larger than 800 GeV. For the graviton interpretation, the boosted analysis is more sensitive than the resolved analysis in the entire studied mass range. For the scalar interpretation, the boosted analysis is more sensitive for m X > 1300 GeV.
The final result is shown in Figure 12. The observed upper limits on the production cross sections range from 5.6 pb for m X = 500 GeV to 0.51 pb for m X = 3000 GeV in the case of a scalar hypothesis. For graviton signals in the same mass interval, they range from 21 pb (18 pb) to 0.28 pb (0.31 pb) in the case of a hypothesis with c = 2.0 (c = 1.0). No boosted analysis was performed for the non-resonant SM signal model.

Conclusion
A search for resonant and non-resonant Higgs boson pair production in the bbWW * decay mode is performed in the bb νqq final state using pp collision data corresponding to an integrated luminosity of 36.1 fb −1 , collected at √ s = 13 TeV by the ATLAS detector at the Large Hadron Collider. No evidence of an excess of events over the background expectation is found. Limits are set on resonant production as a function of the resonance mass for a scalar resonance and for spin-2 gravitons in the mass range 500 to 3000 GeV. An upper limit is set on the cross section of non-resonant pair production σ(pp → HH) · B(HH → bbWW * ) < 2.5 pb at 95% CL corresponding to 300 times the predicted SM cross section. Given the result of this work, in order to bring relevant sensitivity improvement to the HH non-resonant SM searches in this channel at the LHC and at future colliders, more advanced analysis techniques, development of new methods for the normalisation of the tt background, and a more refined estimation of the multijet background, need to be deployed.    [31] T. Sjostrand