Searches for Higgs boson pair production in the $hh\to bb\tau\tau, \gamma\gamma WW*, \gamma\gamma bb, bbbb$ channels with the ATLAS detector

Searches for both resonant and non-resonant Higgs boson pair production are performed in the $hh\to bb\tau\tau, \gamma\gamma WW^*$ final states using 20.3 fb$^{-1}$ of $pp$ collision data at a center-of-mass energy of 8 TeV recorded with the ATLAS detector at the Large Hadron Collider. No evidence of their production is observed and 95\% confidence level upper limits on the production cross sections are set. These results are then combined with the published results of the $hh\to \gamma\gamma bb, bbbb$ analyses. An upper limit of 0.69 (0.47) pb on the non-resonant Standard Model like $hh$ production is observed (expected), corresponding to 70 (48) times of the SM $gg\to hh$ cross section. For production via narrow resonances, cross section limits of $hh$ production from a heavy Higgs boson decay are set as a function of the heavy Higgs boson mass. The observed (expected) limits range from 2.1 (1.1) pb at 260 GeV to 0.011 (0.018) pb at 1000 GeV. These results are interpreted in the context of two simplified scenarios of the Minimal Supersymmetric Standard Model.


Introduction
The Higgs boson discovered at the LHC in 2012 [1, 2] opens a window for testing the scalar sector of the Standard Model (SM) and its possible extensions. Since the discovery, significant progress has been made in measuring its coupling strengths to fermions and vector bosons [3-6] as well as in studying its spin and its charge-conjugate and parity (CP) properties [7,8]. All results are consistent with those expected for the SM Higgs boson (here denoted by h). Within the SM, the existence of the Higgs boson is a consequence of the electroweak symmetry breaking (EWSB). This also predicts self-coupling between Higgs bosons, the measurement of which is crucial in testing the mechanism of EWSB. The self-coupling is one mechanism for Higgs boson pair production as shown in Fig. 1(a). Higgs boson pairs can also be produced through other interactions such as the Higgs-fermion Yukawa interactions ( Fig. 1(b)) in the Standard Model. These processes are collectively referred to as nonresonant production in this paper. Higgs boson pair production at the LHC as a probe of the self-coupling has been extensively studied in the literature [9-13]. One conclusion [14] is that the data collected so far (approximately 25 fb −1 in total) are insensitive to the self-coupling in the SM, because of the expected small signal rates [15][16][17] and large backgrounds. However, it is essential to quantify the sensitivity of the current dataset and to develop tools for future measurements. Moreover, physics beyond the Standard Model (BSM) can potentially enhance the production rate and alter the event kinematics. For example, in the Minimal Supersymmetric Standard Model (MSSM) [18], a heavy CP-even neutral Higgs boson H can decay to a pair of lighter Higgs bosons. Production of H followed by its decay H → hh would lead to a new resonant process of Higgs boson pair production, in contrast to the nonresonant hh production predicted by the SM (Fig. 1). In composite Higgs models such as those discussed in Refs. [19,20], increased production of nonresonant Higgs boson pairs is also expected.
Both the ATLAS and CMS collaborations have searched for nonresonant and/or resonant Higgs boson pair production [21][22][23]. In particular, ATLAS has published the results of searches in the hh → γγbb [21] and hh → bbbb [22] decay channels. 1 In this paper, searches in two additional hh decay final states, bbττ and γγWW * , are reported. For the hh → bbττ analysis, one tau lepton is required to decay to an electron or a muon, collectively referred to as , and the other tau lepton decays to hadrons (τ had ). For hh → γγWW * , the h → WW * → νqq decay signature is considered in this study. The results of these new analyses are combined with the published results of hh → γγbb and hh → bbbb for both nonresonant and resonant production. The resonance mass m H considered in this paper ranges from 260 GeV to 1 Notations indicating particle charges or antiparticles are generally omitted throughout this paper.

Data and Monte Carlo samples
The data used in the searches were recorded in 2012 with the ATLAS detector at the Large Hadron Collider in proton-proton collisions at a center-of-mass energy of 8 TeV and correspond to an integrated luminosity of 20.3 fb −1 . The ATLAS detector is described in detail in Ref. [31]. Only data recorded when all subdetector systems were properly functional are used.
Signal and background MC samples are simulated with various event generators, each interfaced to Pythia v8.175 [32] for parton showers, hadronization and underlying-event simulation. Parton distribution functions (PDFs) CT10 [33] or CTEQ6L1 [34] for the proton are used depending on the generator in question. MSTW2008 [35] and NNPDF [36] PDFs are used to evaluate systematic uncertainties. Table 1 gives a brief overview of the event generators, PDFs and cross sections used for the hh → bbττ and hh → γγWW * analyses. All MC samples are passed through the ATLAS detector simulation program [37] based on GEANT4 [38].
Signal samples for both nonresonant and resonant Higgs boson pair production are generated using the leading-order MadGraph5 v1.5.14 [39] program. The nonresonant production is modeled using the SM DiHiggs model [40,41] while the resonant production is realized using the HeavyScalar model [42], both implemented in MadGraph5. The heavy scalar H is assumed to have a narrow decay width of 10 MeV, much smaller than the experimental resolution. The SM prediction for the nonresonant gg → hh production cross section is 9.9 ± 1.3 fb [17] with m h = 125.4 GeV from the next-to-next-to-leading-order calculation in QCD.
Single SM Higgs boson production is considered as a background. The Powheg r1655 generator [43][44][45] is used to produce gluon fusion (ggF) and vector-boson fusion (VBF) events. This generator calculates QCD corrections up to next-to-leading order (NLO), including finite bottom-and top-quark mass . The Higgs boson transverse momentum (p T ) spectrum of the ggF process is matched to the calculated spectrum at next-to-next-to-leading order (NNLO) and next-to-next-to-leading logarithm (NNLL) [47] in QCD corrections. Events of associated production qq → Vh (here V = W or Z) and qq/gg → tth are produced using the Pythia8 generator [32]. All of these backgrounds are normalized using the state-of-the-art theoretical cross sections (see Table 1 Electrons are reconstructed from energy clusters in the electromagnetic calorimeter matched to tracks in the inner tracker. The calorimeter shower profiles of electron candidates must be consistent with those expected from electromagnetic interactions. Electron candidates are identified using tight and medium criteria [57] for the hh → bbττ and hh → γγWW * analyses, respectively. The selected candidates are required to have transverse momenta 2 p T > 15 GeV and be within the detector fiducial volume of |η| < 2.47 excluding 1.37 < |η| < 1.52, the transition region between the barrel and endcap calorimeters. Muons are identified by matching tracks or segments reconstructed in the muon spectrometer with tracks reconstructed in the inner tracker. They are required to have p T > 10 GeV and |η| < 2.5. Both the electrons and muons must satisfy calorimeter and track isolation requirements. The calorimeter isolation requires that the energy deposited in the calorimeter in a cone of size ∆R ≡ (∆η) 2 + (∆φ) 2 = 0.2 around the lepton (electron or muon), excluding the energy deposited by the lepton itself, is less than 6% (20%) of the p T of the lepton for the hh → bbττ (hh → γγWW * ) analysis. The track isolation is defined similarly: the scalar p T sum of additional tracks originating from the primary vertex with p T > 1 GeV in a cone of size ∆R = 0.4 around the lepton is required to be less than 6% (15%) of the p T of the lepton track for the hh → bbττ (hh → γγWW * ) analysis.
Photons are reconstructed from energy clusters in the electromagnetic calorimeter with their shower profiles consistent with electromagnetic showers. A significant fraction of photons convert into e + e − pairs inside the inner tracker. Thus photon candidates are divided into unconverted and converted categories. Clusters without matching tracks are considered as unconverted photons, while clusters matched to tracks consistent with conversions are considered as converted photons. Photon candidates must fulfill the tight identification criteria [58] and are required to have p T > 25 GeV and |η| < 2.37 (excluding the transition region 1.37 < |η| < 1.52) and must satisfy both calorimeter and track isolation. The calorimeter isolation requires the additional energy in a cone of ∆R = 0.4 around the photon candidate to be less than 6 GeV while the track isolation requires the scalar p T sum of additional tracks in a cone of ∆R = 0.2 around the photon to be less than 2.6 GeV.
Jets are reconstructed using the anti-k t algorithm [59] with a radius parameter of R = 0.4. Their energies are corrected for the average contributions from pileup interactions. Jets are required to have p T > 30 GeV and |η| < 4.5. For the hh → γγWW * analysis, a lower p T requirement of 25 GeV is applied for jets in the central region of |η| < 2.4. To suppress contributions from pileup interactions, jets with p T < 50 GeV and within the acceptance of the inner tracker (|η| < 2.4) must have over 50% of the scalar sum of the p T of their associated tracks contributed by those originating from the primary vertex. Jets containing b-hadrons are identified using a multivariate algorithm (b-tagging) [60]. The algorithm combines information such as the explicit reconstruction of the secondary decay vertices and track impact-parameter significances. The operating point chosen for both hh → bbττ and hh → γγWW * analyses has an efficiency of 80% for the b-quark jets in tt events.
Hadronically decaying τ candidates (τ had ) are reconstructed using clusters in the electromagnetic and hadronic calorimeters [61]. The tau candidates are required to have p T > 20 GeV and |η| < 2.5. The number of tracks with p T > 1 GeV associated with the candidates must be one or three and the total charge determined from these tracks must be ±1. The tau identification uses calorimeter cluster as well as tracking-based variables, combined using a boosted-decision-tree method [61]. Three working points, 2 ATLAS uses a right-hand coordinate system with the interaction point as its origin and the beam line as its z-axis. The x-axis points to the center of the LHC ring and y-axis points upwards. The pseudorapidity η is defined as η = − ln tan(θ/2), where θ is the polar angle measured with respect to the z-axis. The transverse momentum p T is calculated from the momentum p: p T = p sin θ.
labeled loose, medium and tight [61], corresponding to different identification efficiencies are used. Dedicated algorithms that suppress electrons and muons misidentified as τ had candidates are applied as well.
The missing transverse momentum (with magnitude E miss T ) is the negative of the vector sum of the transverse momenta of all photon, electron, muon, τ had , and jet candidates, as well as the p T of all calorimeter clusters not associated with these reconstructed objects, called the soft-term contribution [62]. The hh → bbττ analysis uses the version of the E miss T calculation in the h → ττ analysis [6]. In this calculation, the soft-term contribution is scaled by a vertex fraction, defined as the ratio of the summed scalar p T of all tracks from the primary vertex not matched with the reconstructed objects to the summed scalar p T of all tracks in the event. The hh → γγWW * analysis, on the other hand, uses the E miss T -significance employed by the h → γγ study [56]. It is defined as the ratio of the measured E miss T to its expected resolution estimated using the square root of the scalar sum of the transverse energies of all objects contributing to the E miss T calculation.

Summary of hh → γγbb
The hh → γγbb analysis, published in Ref.
[21], largely follows the ATLAS analysis of the Higgs boson discovery in the h → γγ decay channel [1, 56]. The search is performed in the √ s = 8 TeV dataset corresponding to an integrated luminosity of 20.3 fb −1 . The data were recorded with diphoton triggers that are nearly 100% efficient for events satisfying the photon requirements. Events must contain two isolated photons. The p T for the leading (subleading) photon must be larger than 35% (25%) of the diphoton invariant mass m γγ , which itself must be in the range of 105 < m γγ < 160 GeV. Events must also contain two energetic b-tagged jets; the leading (subleading) jet must have p T > 55 (35) GeV, and the dijet mass must fall within a window 95 < m bb < 135 GeV, as expected from the h → bb decay. A multivariate b-tagging algorithm [60] that is 70% efficient for the b-quark jets in tt events is applied.
Backgrounds for both the resonant and nonresonant analyses are divided into two categories: events containing a single real Higgs boson (with h → γγ), and the continuum background of events not containing a Higgs boson. The former are evaluated purely from simulation, and are small compared to the continuum background, which is evaluated from data in the diphoton mass sidebands (the m γγ range of 105-160 GeV excluding the region of m h ± 5 GeV). In the nonresonant analysis, an unbinned signal-plus-background fit is performed on the observed m γγ distribution, with the background from single Higgs bosons constrained to the expectation from the SM. The continuum background is modeled with an exponential function; the shape of the exponential function is taken from data containing a diphoton and dijet pair where fewer than two jets are b-tagged.
The resonant search proceeds in a similar manner, although it is ultimately a counting experiment, with an additional requirement on the four-object invariant mass m γγbb , calculated with the bb mass constrained to m h . This requirement on m γγbb varies with the resonance mass hypothesis under evaluation, and is defined as the smallest window containing 95% of the signal events based on MC simulation. As in the nonresonant search, the number of background events with real Higgs bosons is estimated from simulation. The continuum background in the m γγ signal window is extrapolated from the diphoton mass sidebands. A resonance with mass between 260 GeV and 500 GeV is considered in the search.
The small number of events (nine) in the diphoton mass sideband leads to large statistical uncertainties (33%) on the dominant continuum background, so that most systematic uncertainties have a small effect on the final result. For the resonant search, however, systematic uncertainties with comparable effect remain. Uncertainties of 0-30% (depending on the resonance mass hypothesis under consideration) are assigned due to the modeling of the m γγbb shape using the data with less than two b-tagged jets. Additional uncertainties of 16-30% arise from the choice of functional form used to parameterize the shape of m γγbb .
In the nonresonant analysis, extrapolating the sidebands into the diphoton mass window of the signal selection leads to a prediction of 1.3 continuum background events. An additional contribution of 0.2 events is predicted from single Higgs boson production. A total of five events are observed, representing an excess of 2.4 standard deviations (σ). A 95% confidence level (CL) upper limit of 2.2 (1.0) pb is observed (expected) for σ(gg → hh), the cross section of nonresonant Higgs boson pair production. For the resonant searches, the observed (expected) upper limits on σ(gg → H) × BR(H → hh) are 2.3 (1.7) pb at m H = 260 GeV and 0.7 (0.7) pb at m H = 500 GeV.

Summary of hh → bbbb
The hh → bbbb analysis [22] benefits from the large branching ratio of h → bb. The analysis employs resolved as well as boosted Higgs boson reconstruction methods. The resolved method attempts to reconstruct and identify separate b-quark jets from the h → bb decay, while the boosted method uses a jet substructure technique to identify the h → bb decay reconstructed as a single jet. The latter is expected if the Higgs boson h has a high momentum. The boosted method is particularly suited to the search for a resonance with mass above approximately 1000 GeV decaying to a pair of SM Higgs bosons. For the combinations presented in this paper, resonances below this mass are considered and the resolved method is used as it is more sensitive.
The analysis with the resolved method searches for two back-to-back and high-momentum bb systems with their masses consistent with m h in a dataset at √ s = 8 TeV corresponding to an integrated luminosity of 19.5 fb −1 for the triggers used. The data were recorded with a combination of multijet triggers using information including the b-quark jet tagging. The trigger is > 99.5% efficient for signal events satisfying the offline selection. Candidate events are required to have at least four b-tagged jets, each with p T > 40 GeV. As in the hh → γγbb analysis, a multivariate b-tagging algorithm [60] with an estimated efficiency of 70% is used to tag jets containing b-hadrons. The four highest-p T b-tagged jets are then used to form two dijet systems, requiring the angular separation ∆R in (η, φ) space between the two jets in each of the two dijet systems to be smaller than 1.5. The transverse momenta of the leading and subleading dijet systems must be greater than 200 GeV and 150 GeV, respectively. These selection criteria, driven partly by the corresponding jet trigger thresholds and partly by the necessity to suppress the backgrounds, lead to significant loss of signal acceptance for lower resonance masses. The resonant search only considers masses above 500 GeV. The leading (m 12 ) and subleading (m 34 ) dijet invariant mass values are required to be consistent with those expected for the hh → bbbb decay, satisfying the requirement: Here m 0 12 (124 GeV) and m 0 34 (115 GeV) are the expected peak values from simulation for the leading and subleading dijet pair, respectively, and σ 12 and σ 34 are the dijet mass resolutions, estimated from the simulation to be 10% of the dijet mass values. More details about the selection can be found in Ref. [22].
After the full selection, more than 90% of the total background in the signal region is estimated to be multijet events, while the rest is mostly tt events. The Z+jets background constitutes less than 1% of the total background and is modeled using MC simulation. The multijet background is modeled using a fully data-driven approach in an independent control sample passing the same selection as the signal except that only one of the two selected dijets is b-tagged. This control sample is corrected for the kinematic differences arising from the additional b-tagging requirements in the signal sample. The tt contribution is taken from MC simulations normalized to data in dedicated control samples.
The dominant sources of systematic uncertainty in the analysis are the b-tagging calibration and the modeling of the multijet background. The degradation in the analysis sensitivity from these uncertainties is below 10%. Other sources of systematic uncertainty include the tt modeling, and the jet energy scale and resolution, which are all at the percent level. A total of 87 events are observed in the data, in good agreement with the SM expectation of 87.0 ± 5.6 events. Good agreement is also observed in the four-jet invariant mass distribution, thus there is no evidence of Higgs boson pair production. For the nonresonant search, both the observed and expected 95% CL upper limit on the cross section σ(pp → hh → bbbb) is 202 fb. For the resonant search, the invariant mass of the four jets is used as the final discriminant from which the upper limit on the potential signal cross section is extracted. The resulting observed (expected) 95% CL upper limit on σ(pp → H → hh → bbbb) ranges from 52 (56) fb, at m H = 500 GeV, to 3.6 (5.8) fb, at m H = 1000 GeV.

hh → bbττ
This section describes the search for Higgs boson pair production in the hh → bbττ decay channel, where only the final state where one tau lepton decays hadronically and the other decays leptonically, bb τ had , is used. The data were recorded with triggers requiring at least one lepton with p T > 24 GeV. These triggers are nearly 100% efficient for events passing the final selection. Candidate bb τ had events are selected by requiring exactly one lepton with p T > 26 GeV, one hadronically decaying tau lepton of the opposite charge with p T > 20 GeV meeting the medium criteria [61], and two or more jets with p T > 30 GeV. In addition, between one and three of the selected jets must be b-tagged using the multivariate b-tagger. The upper bound on the number of b-tagged jets is designed to make this analysis statistically independent of the hh → bbbb analysis summarized in Sec. 5. These criteria are collectively referred to as the "preselection".
The backgrounds from W+jets, Z → ττ, diboson (WW, WZ and ZZ) and top quark (both tt and single top quark) production dominate the surviving sample and their contributions are derived from a mixture of data-driven methods and simulation. The contribution from events with a jet misidentified as a τ had , referred to as the fake τ had background, are estimated using data with a "fake-factor" method. The method estimates contributions from W+jets, multijet, Z+jets and top quark events that pass the event selection due to misidentified τ had candidates. The fake factor is defined as the ratio of the number of τ had candidates identified as medium, to the number satisfying the loose, but not the medium, criteria [61]. The p Tdependent fake factors are measured in data control samples separately for the τ had candidates with one or three tracks and for W+jets, multijet, Z+jets and top quark contributions. The W+jets, multijet, Z+jets and top quark control samples are selected by reversing the m T cut (see below), relaxing the lepton isolation requirement, reversing the dilepton veto or by requiring b-tagged jets, respectively. The fake factors determined from these control samples are consistent within their statistical uncertainties. They are then applied to the signal control sample, i.e., events passing the selection, except that the τ had candidate is required to satisfy the loose, but not the medium, τ had identification, to estimate the fake τ had background. The composition of the sample is determined from a mixture of data-driven methods and MC simulations and it is found that the sample is dominated by the W+jets and multijet events. Details of the method can be found in Ref.
[61]. The method is validated using the same-sign τ had events that are otherwise selected in the same way as the signal candidates.
The Z → ττ background is modeled using selected Z → µµ events in data through embedding [63], where the muon tracks and associated energy depositions in the calorimeters are replaced by the corresponding simulated signatures of the final-state particles of tau decays. In this approach, the kinematics of the produced boson, the hadronic activity of the event (jets and underlying event) as well as pileup are taken from data [6]. Other processes passing the Z → µµ selection, primarily from top quark production, are subtracted from the embedded data sample using simulation. Their contributions are approximately 2% for events with one b-tagged jet and 25% for events with two or more b-tagged jets. The Z → ττ background derived is found to be in a good agreement with that obtained from the MC simulation.
The remaining backgrounds, mostly tt and diboson events with genuine τ had in their decay final states, are estimated using simulation. The small contributions from single SM Higgs boson production and from Z(→ ee/µµ) + jets events (in which one of the electrons or muons is misidentified as τ had ) are also estimated from simulation. The production rates of these processes are normalized to the theoretical cross sections discussed in Sec. 2. For the simulation of the tt process, the top quark p T distribution is corrected for the observed difference between data and simulation [64]. The background from misidentified leptons is found to be negligible.
Figures 2(a) and 2(b) compare the observed ditau (m ττ ) and dijet (m bb ) mass distributions with those expected from background events after the preselection discussed above. The sample is dominated by contributions from top quark production, fake τ had , and Z → ττ backgrounds. Also shown in the figures are the expected signal distributions for a Higgs boson pair production cross section of 20 pb as an illustration. The yield of the nonresonant production is significantly higher than that of the resonant production for the same cross section, largely due to the harder p T spectrum of the Higgs boson h of the nonresonant production. The ditau invariant mass is reconstructed from the electron or muon, τ had , and E miss T using a method known as the missing mass calculator (MMC) [65]. The MMC solves an underconstrained system of equations with solutions weighted by E miss T resolution and the tau-lepton decay topologies. It returns the most probable value of the ditau mass, assuming that the observed lepton, τ had and E miss T stem from a ττ resonance. The dijet mass is calculated from the two leading b-tagged jets, or using also the highest-p T untagged jet if only one jet is b-tagged.
Additional topological requirements are applied to reduce the background. As shown in Fig. 2(c), the signal events tend to have small values of the transverse mass m ν T calculated from the lepton and E miss T system. Consequently, a requirement of m ν T < 60 GeV is applied, which reduces the background significantly with only a small loss of the signal efficiency. In addition, the angular separation in the transverse plane between the E miss T and τ had is required to be larger than one radian to reduce the fake τ had background.
Background events from tt → WWbb → ντνbb decay have an identical visible final state to the signal if the tau lepton decays hadronically. For signal h → ττ → τ had events, however, the p T of the lepton tends to be softer than that of the τ had due to the presence of more neutrinos in the τ → decays. Thus the p T of the electron or muon is required to satisfy p T ( ) < p T (τ had ) + 20 GeV. The tt events of the tt → WW bb → ν qq bb final state with a misidentified τ had candidate remain a large background. To reduce its contribution, a W boson candidate is reconstructed from the τ had candidate and its closest system. The top quark background includes contributions from both tt and the single top-quark production. The background category labeled "Others" comprises diboson and Z → ee/µµ contributions. Contributions from single SM Higgs boson production are included in the background estimates, but are too small to be visible on these distributions. As illustrations, expected signal distributions for a Higgs boson pair production cross section of 20 pb are overlaid for both nonresonant and resonant Higgs boson pair production. A mass of m H = 300 GeV is assumed for the resonant production. The last bin in all distributions represents overflows. The gray hatched bands represent the uncertainties on the total background contributions. These uncertainties are largely correlated from bin to bin. untagged jet and its mass m τ j is calculated. The W candidate is then paired with a b-tagged jet to form the top quark candidate with a reconstructed mass m τ jb . The pairing is chosen to minimize the mass sum m b + m τb for events with two or more b-tagged jets. If only one jet is b-tagged, one of the b-jets in the mass sum is replaced by the highest-p T untagged jet. An elliptical requirement in the form of a χ 2 in the (m τ j , m τ jb ) plane: ∆m W cos θ − ∆m t sin θ 28 GeV 2 + ∆m W sin θ + ∆m t cos θ 18 GeV 2 > 1 is applied to reject events with (m τ j , m τ jb ) consistent with the hypothesis (m W , m t ), the masses of the W boson and the top quark. Here ∆m W = m τ j − m W , ∆m t = m τ jb − m t and θ is a rotation angle given by tan θ = m t /m W to take into account the average correlation between m τ j and m τ jb .
Finally, the remaining events must have 90 < m bb < 160 GeV, consistent with the expectation for the h → bb decay. For the nonresonant search, m ττ is used as the final discriminant to extract the signal, and its distribution is shown in Fig. 3(a). The selection efficiency for the gg → hh → bbττ signal is 0.57%. For the resonant search, the MMC mass is required to be in the range of 100 < m ττ < 150 GeV. The mass resolutions of m bb and m ττ are comparable for the signal, but the m bb distribution has a longer tail. The resonance mass m bbττ reconstructed from the dijet and ditau system is used as the discriminant. To improve the mass resolution of the heavy resonances, scale factors of m h /m bb and m h /m ττ are applied respectively to the four-momenta of the bb and ττ systems, where m h is set to the value of 125 GeV used in the simulation. The resolution of m bbττ is found from simulation studies to vary from 2.4% at m H = 260 GeV to 4.8% at 1000 GeV. The improvement in the resolution from the rescaling is largest at low mass and varies from approximately a factor of three at 260 GeV to about 30% at 1000 GeV. The reconstructed m bbττ distribution for events passing all the selections is shown in Fig. 3(b). The efficiency for the gg → H → hh → bbττ signal varies from 0.20% at 260 GeV to 1.5% at 1000 GeV. These efficiencies include branching ratios of the tau lepton decays, but not those of heavy or light Higgs bosons.
To take advantage of different signal-to-background ratios in different kinematic regions, the selected events are divided into four categories based on the ditau transverse momentum p ττ T (less than or greater than 100 GeV) and the number of b-tagged jets (n b = 1 or ≥ 2) for both the nonresonant and resonant searches. The numbers of events expected from background processes and observed in the data passing the resonant hh → bbττ selection are summarized in Table 2 for each of the four categories. The expected number of events from the production of a Higgs boson with m H = 300 GeV and a cross section of σ(gg → H) × BR(H → hh) = 1 pb for each category is also shown for comparison.
Systematic uncertainties from the trigger, luminosity, object identification, background estimate as well as Monte Carlo modeling of signal and background processes are taken into account in the background estimates and the calculation of signal yields. The impact of these systematic uncertainties varies for different background components and event categories. For the most sensitive n b ≥ 2 categories, the main background contributions are from top quark, fake τ had , and Z → ττ. The jet energy scale and resolution is the largest uncertainty for the top-quark contribution, ranging between 10% and 19% for the nonresonant and resonant searches. The leading source of systematic uncertainty for the fake τ had background is the "fake-factor" determination, due to the uncertainties of the sample composition. Varying the composition of W+jets, Z+jets, top quark and multijet events in the control samples by ±50% leads to a change in the estimated fake τ had background by 9.5%. The most important source of systematic uncertainty for the Z → ττ background is the tt subtraction from the Z → µµ sample used for the embedding, due to the uncertainty on the tt normalization. Its effect ranges from 8% to 15%. The overall systematic uncertainties on the total background contributions to the high (low) p ττ T category of n b ≥ 2 are 12% (9%) for the nonresonant search and 14% (14%) for the resonant search. The largest contributions are Table 2: The numbers of events predicted from background processes and observed in the data passing the final selection of the resonant search for the four categories. The top quark background includes contributions from both tt and the single top-quark production. The "others" background comprises diboson and Z → ee/µµ contributions. The numbers of events expected from the production of a m H = 300 GeV Higgs boson with a cross section of σ(gg → H) × BR(H → hh) = 1 pb are also shown as illustrations. The uncertainties shown are the total uncertainties, combining statistical and systematic components.  Contributions from single SM Higgs boson production are included in the background estimates, but are too small to be visible on these distributions. As illustrations, the expected signal distributions assume a cross section of 10 pb for Higgs boson pair production for both the nonresonant and resonant searches. In (b), a resonance mass of 300 GeV is assumed. The gray hatched bands represent the uncertainties on the total backgrounds. These uncertainties are largely correlated from bin to bin. from jet and tau energy scales and b-tagging. The modeling of top quark production is also an important systematic uncertainty for the category with two or more b-tagged jets and high p ττ T . The uncertainties on the signal acceptances are estimated from experimental as well as theoretical sources. The total experimental systematic uncertainties vary between 12% and 24% for the categories with two or more b-tagged jets, and are dominated by the jet and tau energy scales and b-tagging. Theoretical uncertainties arise from the choice of parton distribution functions, the renormalization and factorization scales as well as the value of strong coupling constant α s used to generate the signal events. Uncertainties of 3%, 1% and 3% from the three sources, respectively, are assigned to all signal acceptances.
For the nonresonant search, the observed ditau mass distribution agrees well with that of the estimated background events as shown in Fig. 3(a). For the resonant search, a small deficit with a local significance of approximately 2σ is observed in the data relative to the background expectation at m bbττ ∼ 300 GeV as is shown in Fig. 3(b). No evidence of Higgs boson pair production is present in the data. The resulting upper limits on Higgs boson pair production from these searches are described in Sec. 9. 7 hh → γγWW * This section describes the search for Higgs boson pair production in the hh → γγWW * decay channel, where one Higgs boson decays to a pair of photons and the other decays to a pair of W bosons. The h → γγ decay is well suited for tagging the Higgs boson. The small Higgs boson width together with the excellent detector resolution for the diphoton mass strongly suppresses background contributions. Moreover, the h → WW * decay has the largest branching ratio after h → bb. To reduce multijet backgrounds, one of the W bosons is required to decay to an electron or a muon (either directly or through a tau lepton) whereas the other is required to decay hadronically, leading to the γγ νqq final state.
The data used in this analysis were recorded with diphoton triggers with an efficiency close to 100% for diphoton events passing the final offline selection. The diphoton selection follows closely that of the ATLAS measurement of the h → γγ production rate [56] and that of the hh → γγbb analysis [21]. Events are required to have two or more identified photons with the leading and subleading photon candidates having p T /m γγ > 0.35 and 0.25, respectively, where m γγ is the invariant mass of the two selected photons. Only events with m γγ in the range of 105 < m γγ < 160 GeV are considered.
Additional requirements are applied to identify the h → WW * → νqq decay signature. Events must have two or more jets, and exactly one lepton, satisfying the identification criteria described in Sec. 3. To reduce multijet backgrounds, the events are required to have E miss T with significance greater than one. Events with any b-tagged jet are vetoed to reduce contributions from top quark production.
A total of 13 events pass the above selection. The final hh → γγWW * candidates are selected by requiring the diphoton mass m γγ to be within a ±2σ window of the Higgs boson mass in h → γγ where σ is taken to be 1.7 GeV. Due to the small number of events, both nonresonant and resonant searches proceed as counting experiments. The selection efficiency for the hh → γγWW * signal of SM nonresonant Higgs boson pair production is estimated using simulation to be 2.9%. For the resonant production, the corresponding efficiency varies from 1.7% at 260 GeV to 3.3% at 500 GeV. These efficiencies include the branching ratios of the W boson decays, but not those of the Higgs boson decays.
The background contributions considered are single SM Higgs boson production (gluon fusion, vectorboson fusion and associated production of Wh, Zh and tth) and continuum backgrounds in the m γγ spectrum. Events from single Higgs boson production can mimic the hh → γγWW * signal if, for example, the Higgs boson decays to two photons and the rest of the event satisfies the h → WW * → νqq identification. These events would exhibit a diphoton mass peak at m h . As in the hh → bbττ analysis, their contributions are estimated from simulation using the SM cross sections [27]. The systematic uncertainty on the total yield of these backgrounds is estimated to be 29%, dominated by the modeling of jet production (27%). The total number of events expected from single SM Higgs production is therefore 0.25 ± 0.07 with contributions of 0.14, 0.08 and 0.025 events from Wh, tth and Zh processes, respectively. Contributions from gluon and vector-boson fusion processes are negligible. The background that is nonresonant in the γγ mass spectrum is measured using the continuum background in the m γγ spectrum. The major source of these backgrounds is Wγγ + jets events with a W → ν decay. These events are expected to have a diphoton mass distribution with no resonant structure at m h and their contribution (N est SR ) in the signal region, m γγ ∈ m h ± 2σ, is estimated from the m γγ sidebands in the data: Here N Data SB is the number of events in the data sidebands, defined as the mass region 105 < m γγ < 160 GeV excluding the signal region. The quantity f SB is the fraction of background events in 105 < m γγ < 160 GeV falling into the signal mass window, and can in principle be determined from a fit of the observed m γγ distribution to an ansatz function. However, the small number of events after the final selection makes such a fit unsuitable. Instead, f SB is determined in a data control sample, selected as the signal sample without the lepton and E miss T requirements. Figure 4(a) shows the m γγ distribution of events in the control sample. For the fit, an exponential function is used to model the sidebands and a wider region of m h ± 5 GeV is excluded to minimize potential signal contamination in the sidebands. The fit yields a value of f SB = 0.1348 ± 0.0001. Varying the fit range of the sidebands leads to negligible changes. Different fit functions, such as a second-order polynomial or an exponentiated second-order polynomial, lead to a difference of 1.4% in f SB . To study the sample dependence of f SB , the fit is repeated for the control sample without the jet and E miss T requirements and a difference of only 2% is observed. Simulation studies show that the continuum background is dominated by W( ν)γγ + jets production. The γγ ν + jets events generated using MadGraph reproduce well the observed m γγ distribution. The potential difference between γγ +jets and γγ ν+jets samples is studied using simulation. A difference below 1% is observed. Taking all these differences as systematic uncertainties, the fraction of background events in the signal mass window is f SB = 0.135 ± 0.004. With 9 (N Data SB ) events observed in the data sidebands, it leads to N est SR = 1.40 ± 0.47 events from the continuum background. Figure 4(a) also shows the contribution expected from single SM Higgs boson production. The data prefer a larger cross section than the SM prediction for single Higgs boson production, consistent with the measurement reported in Ref. [66].
The uncertainties on the signal acceptances are estimated following the same procedure as the hh → bbττ analysis. The total experimental uncertainty is found to vary between 4% and 7% for different signal samples under consideration, dominated by the contribution from the jet energy scale. The theoretical uncertainties from PDFs, the renormalization and factorization scales, and the strong coupling constant are 3%, 1%, and 3%, respectively, the same as for the hh → bbττ analysis.
The m γγ distribution of the selected events in the data is shown in Fig. 4(b). In total, 13 events are found with 105 < m γγ < 160 GeV. Among them, 4 events are in the signal mass window of m h ± 2σ compared with 1.65 ± 0.47 events expected from single SM Higgs boson production and continuum background processes. The p-value of the background-only hypothesis is 3.8%, corresponding to 1.8 standard deviations.
Assuming a cross section of 1 pb (σ(gg → hh) or σ(gg → H) × BR(H → hh)) for Higgs boson pair production, the expected number of signal events is 0.64 ± 0.05 for the nonresonant production. For the resonant production, the corresponding numbers of events are 0.47 ± 0.05 and 0.72 ± 0.06 for a resonance mass of 300 GeV and 500 GeV, respectively. The implications of the search for Higgs boson pair production are discussed in Sec. 9.

Combination procedure
The statistical analysis of the searches is based on the framework described in Refs. [67][68][69][70]. Profilelikelihood-ratio test statistics are used to measure the compatibility of the background-only hypothesis with the observed data and to test the hypothesis of Higgs boson pair production with its cross section as the parameter of interest. Additional nuisance parameters are included to take into account systematic uncertainties and their correlations. The likelihood is the product of terms of the Poisson probability constructed from the final discriminant and of nuisance parameter constraints with either Gaussian, lognormal, or Poisson distributions. Upper limits on the Higgs boson pair production cross section are derived using the CL s method [71]. For the combinations, systematic uncertainties that affect two or more analyses (such as those of luminosity, jet energy scale and resolutions, b-tagging, etc.) are modeled with common nuisance parameters.
For the hh → bbττ analysis, Poisson probability terms are calculated for the four categories separately from the mass distributions of the ditau system for the nonresonant search ( Fig. 3(a)) and of the bbττ system for the resonant search ( Fig. 3(b)). The m bbττ distributions of the resonant search are rebinned to ensure a sufficient number of events for the background prediction in each bin, in particular a single bin is used for m bbττ 400 GeV for each category. For the hh → γγWW * analysis, event yields are used to calculate Poisson probabilities without exploiting shape information. The hh → γγbb and hh → bbbb analyses are published separately in Refs. [21,22]. However, the results are quoted at slightly different values of the Higgs boson mass m h and, therefore, have been updated using a common mass value of m h = 125.4 GeV [24] for the combinations. The decay branching ratios of the Higgs boson h and their uncertainties used in the combinations are taken from Ref. [27]. Table 3 is a summary of the number of categories and final discriminants used for each analysis. The four individual analyses are sensitive to different kinematic regions of the hh production and decays. The combination is performed assuming that the relative contributions of these regions to the total cross section are modeled by the MadGraph5 [39] program used to simulate the hh production.

Results
In this section, the limits on the nonresonant and resonant searches are derived. The results of the hh → bbττ and hh → γγWW * analyses are first determined and then combined with previously published results of the hh → γγbb and hh → bbbb analyses. The impact of the leading systematic uncertainties is also discussed.
The observed and expected upper limits at 95% CL on the cross section of nonresonant production of a Higgs boson pair are shown in Table 4. These limits are to be compared with the SM prediction of 9.9 ± 1.3 fb [17] for gg → hh production with m h = 125.4 GeV. Only the gluon fusion production process is considered. The observed (expected) cross-section limits are 1.6 (1.3) pb and 11.4 (6.7) pb from the hh → bbττ and hh → γγWW * analyses, respectively. Also shown in the table are the cross-section limits relative to the SM expectation. The results are combined with those of the hh → γγbb and hh → bbbb analyses. The p-value of compatibility of the combination with the SM hypothesis is 4.4%, equivalent to 1.7 standard deviations. The low p-value is a result of the excess of events observed in the hh → γγbb analysis. The combined observed (expected) upper limit on σ(gg → hh) is 0.69 (0.47) pb, corresponding to 70 (48) times the cross section predicted by the SM. The hh → bbbb analysis has the best expected sensitivity followed by the hh → γγbb analysis. The observed combined limit is slightly weaker than that of the hh → bbbb analysis, largely due to the aforementioned excess. The impact of systematic uncertainties on the cross-section limits is studied using the signal-strength parameter µ, defined as the ratio of the extracted to the assumed signal cross section (times branching ratio BR(H → hh) for the resonant search). The resulting shifts in µ depend on the actual signal-strength value. For illustration, they are evaluated using a cross section of 1 pb for gg → (H →)hh, comparable to the limits set. The effects of the most important uncertainty sources are shown in Table 5. The leading contributions are from the background modeling, b-tagging, the h decay branching ratios, jet and E miss T measurements. The large impact of the b-tagging systematic uncertainty reflects the relatively large weight of the hh → bbbb analysis in the combination. For the hh → bbττ analysis alone, the three leading systematic sources are the background estimates, jet and E miss T measurements, and lepton and τ had identifications. For the hh → γγWW * analysis, they are the background estimates, jet and E miss T measurements and theoretical uncertainties of the decay branching ratios of the Higgs boson h.
For the resonant production, limits are set on the cross section of gg → H production of the heavy Higgs boson times its branching ratio BR(H → hh) as a function of the heavy Higgs boson mass m H . The observed (expected) limits of the hh → bbττ and hh → γγWW * analyses are illustrated in Fig. 5 and listed Table 5: The impact of the leading systematic uncertainties on the signal-strength parameter µ of a hypothesized signal for both the nonresonant and resonant (m H = 300, 600 GeV) searches. For the signal hypothesis, a Higgs boson pair production cross section (σ(gg → hh) or σ(gg → H) × BR(H → hh)) of 1 pb is assumed. in Table 6 (along with results from the hh → γγbb and hh → bbbb analyses). The m H search ranges are 260-1000 GeV for hh → bbττ and 260-500 GeV for hh → γγWW * . For the hh → bbττ analysis, the observed limit around m H ∼ 300 GeV is considerably lower than the expectation, reflecting the deficit in the observed m bbττ distribution. At high mass, the limits are correlated since a single bin is used for m bbττ 400 GeV. The decrease in the limit as m H increases is a direct consequence of increasing selection efficiency for the signal. This is also true for the hh → γγWW * analysis as the event selection is independent of m H .  The hh → γγbb and hh → bbbb analyses are published separately and the mass range covered by the two analyses are 260-500 GeV and 500-1500 GeV, respectively. The results of these four analyses, summarized in Table 6, are combined for the mass range 260-1000 GeV assuming the SM values of the h decay branching ratios. To reflect the better mass resolutions of the hh → bbbb and hh → γγbb analyses, the combination is performed with smaller mass steps than those of the hh → bbττ and hh → γγWW * analyses. The most significant excess in the combined results is at a resonance mass of 300 GeV with a local significance of 2.5σ, largely due to the 3.0σ excess observed in the hh → γγbb analysis [21]. The upper limit on σ(gg → H) × BR(H → hh) varies from 2.1 pb at 260 GeV to 0.011 pb at 1000 GeV. These limits are shown in Fig. 6 as a function of m H . For the low-mass region of 260-500 GeV, both the hh → γγbb and hh → bbττ analyses contribute significantly to the combined sensitivities. Above 500 GeV, the sensitivity is dominated by the hh → bbbb analysis. Table 5 shows the impact of the leading systematic uncertainties for a heavy Higgs boson mass of 300 GeV and 600 GeV. As in the nonresonant search, the systematic uncertainties with the largest impact on the sensitivity are from the uncertainties on the background modeling, b-tagging, jet and E miss T measurements, and the h decay branching ratios. These limits are directly applicable to models such as those of Refs. [72][73][74][75][76][77] in which the Higgs boson h has the same branching ratios as the SM Higgs boson. uncertainty ranges of the expected combined limits. The improvement above m H = 500 GeV is due to the sensitivity of the hh → bbbb analysis. The more finely spaced mass points of the combination reflect the better mass resolutions of the hh → γγbb and hh → bbbb analyses than those of the hh → bbττ and hh → γγWW * analyses.

Interpretation
The upper cross-section limits of the resonant search are interpreted in two MSSM scenarios, one referred to as the hMSSM [28,29] and the other as the low-tb-high [30]. In the interpretation, the CP-even light and heavy Higgs bosons of the MSSM are assumed to be the Higgs bosons h and H of the search, respectively. The natural width of the heavy Higgs boson H where limits are set in these scenarios is sufficiently smaller than the experimental resolution, which is at best 1.5%, that its effect can be neglected.
In the hMSSM scenario, the mass of the light CP-even Higgs boson is fixed to 125 GeV in the whole parameter space. This is achieved by implicitly allowing the supersymmetry-breaking scale m S to be very large, which is especially true in the low tan β region where m S 1 TeV, and making assumptions about the CP-even Higgs boson mass matrix and its radiative corrections, as well as the Higgs boson coupling dependence on the MSSM parameters. Here tan β is the ratio of the vacuum expectation values of the two doublet Higgs fields. The "low-tb-high" MSSM scenario follows a similar approach, differing in that explicit choices are made for the supersymmetry-breaking parameters [30]. The mass of the light Higgs boson is not fixed in this scenario, but is approximately 125 GeV in most of the parameter space. The m h value grows gradually from 122 GeV at m A ∼ 220 GeV to 125 GeV as m A approaching infinity. Higgs boson production cross sections through the gluon-fusion process are calculated with SusHi 1.4.1 [78][79][80] for both scenarios. Higgs boson decay branching ratios are calculated with HDECAY 6.42 [81] following the prescription of Ref. [29] for the hMSSM scenario and with FeynHiggs 2.10.0 [82][83][84] for the low-tbhigh scenario.  Table 6. The upper limits, as functions of m H , are recomputed; the hh decay fractions for each final state are fixed to their smallest value found in 1 < tan β < 2, the range of the expected sensitivity. This approach yields conservative limits, but simplifies the computation as the limit calculation does not have to be repeated at each tan β value. The results are used to set exclusions in the (tan β, m A ) plane as shown in Fig. 7. The analysis is sensitive to the region of low tan β and m A values in the range ∼ 200-350 GeV. For m A 200 GeV, m H is typically below the 2m h threshold of the H → hh decay, whereas above 350 GeV, the H → hh decay is suppressed because of the dominance of the H → tt decay. The observed exclusion region in the (tan β, m A ) plane is smaller than the expectation, reflecting the small excess observed in the data.

Summary
This paper summarizes the search for both nonresonant and resonant Higgs boson pair production in proton-proton collisions from approximately 20 fb −1 of data at a center-of-mass energy of 8 TeV recorded by the ATLAS detector at the LHC. The search is performed in hh → bbττ and γγWW * final states. No significant excess is observed in the data beyond the background expectation. Upper limits on the hh production cross section are derived. Combining with the hh → γγbb, bbbb searches, a 95% CL upper  [51] N. Kidonakis, Next-to-next-to-leading-order collinear and soft gluon corrections for t-channel single top quark production, Phys. Rev. D83 (2011) 091503, arXiv:1103.2792 [hep-ph].
[64] ATLAS Collaboration, Measurements of normalized differential cross sections for tt production in pp collisions at √ s = 7 TeV using the ATLAS detector, Phys.