A search for Higgs boson pair production in the bb tau tau final state in proton-proton collisions at sqrt(s) = 8 TeV

Results are presented from a search for production of Higgs boson pairs (HH) where one boson decays to a pair of b quarks and the other to a tau lepton pair. This work is based on proton-proton collision data collected by the CMS experiment at sqrt(s) = 8 TeV, corresponding to an integrated luminosity of 18.3 inverse femtobarns. Resonant and nonresonant modes of HH production have been probed and no significant excess relative to the background-only hypotheses has been found in either mode. Upper limits on cross sections of the two HH production modes have been set. The results have been combined with previously published searches at sqrt(s) = 8 TeV, in decay modes to two photons and two b quarks, as well as to four b quarks, which also show no evidence for a signal. Limits from the combination have been set on resonant HH production by an unknown particle X in the mass range m[X] = 300 GeV to m[X] = 1000 GeV. For resonant production of spin 0 (spin 2) particles, the observed 95% CL upper limit is 1.13 pb (1.09 pb) at m[X] = 300 GeV and to 21 fb (18 fb) at m[X] = 1000 GeV. For nonresonant HH production, a limit of 43 times the rate predicted by the standard model has been set.


Introduction
The discovery of a standard model (SM)-like Higgs (H) boson [1,2] motivates further investigation of the nature of electroweak symmetry breaking. In particular, the measurement of the Higgs self-coupling can provide valuable information about the details of the mechanism by which the electroweak symmetry is broken.
The measurement of the H pair (HH) production rate allows us to probe the trilinear H selfcoupling. The leading-order (LO) Feynman diagrams for SM HH production are shown in Fig. 1. The amplitude of the triangle diagram depends on the trilinear H self-coupling. Interference of the box diagram with the triangle diagram reduces the SM cross section to a value of about 10 fb at a center-of-mass energy of √ s = 8 TeV [3]. A deviation of the trilinear H self-coupling from the SM value may enhance the HH production rate significantly. The composite Higgs models discussed in Refs. [4,5] predict such an enhancement in which the mass distribution of the H pair is expected to be broad. We refer to this case as non-resonant HH production.
Alternatively, the HH production rate could be enhanced if a unknown heavy particle X decays into a pair of H's. The LO process for this case is shown in Fig. 2. We refer to this case as resonant HH production. Several models beyond the SM give rise to such decays, in particular, two-Higgs-doublet models [6,7], composite Higgs boson models [4,8], Higgs portal models [9,10], and models involving warped extra dimensions (WED) [11]. The present search is performed in the context of the latter models in which the heavy resonance X can either be a radion with spin 0 [12][13][14][15] or a Kaluza-Klein (KK) excitation of the graviton with spin 2 [16,17]. The benchmark points for both models can be expressed in terms of the dimensionless quantity k/M Pl and the mass scale Λ R = √ 6e −kl M Pl , where k is the exponential warp factor for the extra dimension, l is the size of the extra dimension, and M Pl is the reduced Planck mass which, is defined by M Pl / √ 8π, where M Pl is the Planck mass. The mass scale Λ R is interpreted as the ultraviolet cutoff of the model [18,19]. In this paper we assume that the SM particles within such a theory follow the characteristics of the SM gauge group and that the right-handed top quark is localized on the TeV brane, referred to as the elementary top hypothesis [20]. A possible mixing between the radion and the H (r/H mixing) [21] is neglected, since precision electroweak studies show that the mixing is most likely to be small [22].
Searches for HH production have been performed previously by the CMS Collaboration at the CERN LHC [23][24][25][26][27] in multi-lepton, multi-lepton+γγ, bbττ, γγbb, and bbbb final states. In this paper we present the results for HH production when one of the H's decays to two bottom quarks, and the other decays to two τ leptons, where the τ leptons decay to hadrons and a ν τ (τ h ). This decay channel is important because of its large branching fraction. A previous search in this channel was performed in the mass range of m X = 260-350 GeV [24]. The present work extends that search to a larger range of resonance mass and to the case of non-resonant HH production. The sensitivity of the analysis is enhanced by reconstructing the full four-vector of the H that decays into τ leptons with a likelihood based algorithm and identifying hadronic τ decays with a multivariate algorithm. We combine the results of the search in the bbττ decay channel with those from searches in the γγbb and and bbbb final states in order to increase the sensitivity to potential signals.
The ATLAS Collaboration has searched for resonant as well as non-resonant HH production in the bbττ, γγWW * , γγbb, and bbbb decay channels [28][29][30]. Their observed (expected) limit on non-resonant HH production, obtained by combining all channels, corresponds to 70 (48) times the SM production rate. The observed (expected) limit on non-resonant HH production obtained from the bbτ h τ h channel alone is 160 (130) times the rate expected in the SM. In case of resonant HH production, the ATLAS Collaboration has set a combined observed (expected) limit on the production rate (σ(pp → X) B(X → HH)) that ranges from 2.1 pb (1.1 pb) at m X = 300 GeV to 11 fb (18 fb) at m X = 1000 GeV. The observed (expected) limit set in the bbτ h τ h channel alone ranges from 1.7 pb (3.1 pb) at m X = 300 GeV to 0.46 pb (0.28 pb) at m X = 1000 GeV.

Experimental setup, data, and Monte Carlo events
This Section briefly describes the CMS detector, emphasizing the tracking detector which plays an important role in this analysis. Details of the experimental data set and the Monte Carlo (MC) simulated event samples for signal events as well as various background processes that are relevant to HH production and decay are also given here.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the superconducting volume are a silicon tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. In the tracker the inner 3 (2) layers in the barrel (endcap) region consist of pixel detectors. The outer 10 (12) layers in the barrel (endcap) region are made of strip detectors. The tracker provides a resolution of ∼0.5% for the measurement of transverse momentum (p T ) of tracks which is important for the search described here. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [31]. The CMS trigger system is composed of two levels [32]. The first level, composed of custom hardware processors, reduces the event rate from 40 MHz to 0.1 MHz. At the next stage, the high-level software-based trigger, implemented in a farm of about 10 000 commercial processor cores, reduces the rate further to less than 1 kHz.

Data and simulated samples
This search is based on proton-proton (pp) collision data corresponding to an integrated luminosity of 18.3 fb −1 recorded at √ s = 8 TeV in 2012. On average, 21 inelastic pp interactions per LHC bunch crossing occurred during this period [33]. One of the interactions is selected as the primary interaction and the rest are called "pileup". Signal samples for both resonant and non-resonant HH production are generated with MADGRAPH 5.1 [34]. For resonant HH production, simulated samples are generated for spin 0 (radion) and spin 2 (graviton) hypotheses for the X resonance at masses m X = 300, 500, 700, and 1000 GeV. Shape templates for the mass parameter of the HH system used in the signal extraction procedure described in Section 8 are produced for intermediate mass points using a horizontal template morphing technique [35] in steps of 50 GeV between 300 and 700 GeV mass points and in steps of 100 GeV between 700 and 1000 GeV mass points. The efficiency and the acceptance are interpolated linearly between the mass points.
The background contribution from multijet events is estimated from data, as described in Section 6.1. Background events arising from Z/γ * → ( = e, µ), W+jets, tt, single top quark, and di-boson (WW, WZ, ZZ) production are modeled using MC samples. Among these backgrounds Z/γ * → , W+jets, tt, and di-boson samples are generated with MADGRAPH 5.1, while the single top quark samples are modeled with POWHEG 1.0 [36].
The Z/γ * → and the W+jets backgrounds are generated in bins of generator-level parton multiplicity in order to enhance the event statistics in regions of high signal purity. These samples are normalized to their respective NNLO cross sections [37]. The tt sample is normalized to the top quark pair production cross section measured by CMS [38] multiplied by a correction factor obtained from a tt enriched control region in data. Furthermore, a kinematic reweighting is applied to simulated tt events [39,40] to match the top quark p T distribution observed in data. The single top quark and the di-boson events are normalized to their respective next-toleading order (NLO) cross sections [41].
Production of events with a single H in the SM scenario is simulated using POWHEG 1.0. The production processes considered are gluon-gluon fusion (ggH), vector boson fusion (qqH), associated production of the H with W and Z bosons (VH), bb or tt pairs. These samples are produced for a H of mass m H = 125 GeV and are normalized to the corresponding cross section given in Ref. [42]. The H decays that have been taken into account in this analysis are H → bb for VH production, H → ττ for VH and ggH production, and both H → bb and H → ττ for qqH production.
Parton shower and hadronization processes are modeled using PYTHIA 6.4. Taus are decayed by TAUOLA 27.121.5 [43]. Pileup interactions represented by minimum bias events generated with PYTHIA 6.4 [44] are added to all simulated samples according to the pileup profile observed in data during the 2012 data-taking period. The generated events are passed through a GEANT4 [45] based simulation of the CMS detector and are reconstructed using the same version of the CMS software as that for data.
A special technique, referred to as embedding, is used to model the background arising from Z/γ * → ττ production. Embedded samples are produced by selecting Z/γ * → µµ events in data and replacing the reconstructed muons by generator-level τ leptons with the same fourvectors as that of the muons [46]. The τ lepton decays are simulated using TAUOLA 27.121.5 and their polarization effects are modeled with TAUSPINNER (TAUOLA++ 1.1.4) [47]. The visible decay products of the τ are reconstructed with the particle-flow (PF) algorithm (cf. Section 3), and then added to the remaining particles of the Z/γ * → µµ event, after removing the two muons. Finally, the τ h candidates, the jets, and the missing transverse momentum vector p miss T , which is defined as the negative vectorial sum of the p T of all reconstructed particles, are reconstructed, and the event is analyzed as if it were data.
The sample of Z/γ * → µµ events that is used as input for the production of Z/γ * → ττ embedded samples contains contributions from the background tt While the overall level of this contribution is small (∼0.1% of the Z/γ * → ττ embedded sample), the contamination of the embedded sample with these events becomes relevant for events selected with one or more jets originating from b quarks. The tt contamination is corrected using simulated tt events that are fed through the same embedding procedure as described above.

Physics object reconstruction and identification
This section describes the methods employed to identify various particles used in this analysis. The PF algorithm is used to reconstruct and identify individual particles (referred to as candidates), such as electrons, muons, photons, charged and neutral hadrons with an optimized combination of information from various elements of the CMS detector [48]. The resulting candidates are used to reconstruct jets, hadronic τ decays, and p miss T . It is required that all candidates in an event originate from a common interaction point, the primary vertex. The sum of p 2 T of all tracks associated with each interaction vertex is computed and the one with the largest value is selected as the primary vertex.

Jets and p miss
T Jets within |η| < 4.7 are built using the anti-k T algorithm [49] implemented in the FASTJET package [50], with distance parameter of 0.5, using PF candidates as input. Misreconstructed jets, mainly arising from calorimeter noise, are rejected by requiring the jets to pass a set of loose identification criteria [51]. Jets originating from pileup interactions are suppressed by an identification discriminant [52] based on multivariate (MVA) techniques. Corrections based on the median energy density per event [53,54] as computed by the FASTJET algorithm, are applied to the jet energy in order to correct for other pileup effects. The energy of reconstructed jets is calibrated as a function of p T and η of the jet [55]. Jets of |η| < 2.4 and p T > 20 GeV are tagged as b quark jets if they are selected by an MVA based algorithm which uses lifetime information of b quarks ("Combined Secondary Vertex", CSV, algorithm). The b tagging efficiency and mistag (misidentification of jets without b quarks as b quark jets) rates for this search are 70% and 1.5% (10%) for light (charm) quarks respectively [56].
The magnitude and direction of the p miss T vector are reconstructed using an MVA based algorithm [33] which uses the fact that pileup predominantly produces low-p T jets and 'unclustered energy' (hadrons not within jets), while isolated leptons and high-p T jets are almost exclusively produced by the hard-scatter interaction, even in high-pileup conditions. In addition, the algorithm provides event-by-event estimate of the p miss T resolution.

Lepton identification
Electrons and muons are used in this analysis solely for the purpose of vetoing events, as described in Section 4. A description of the electron and the muon identification criteria and the computation of their isolation from other particles is given in Refs. [57,58].
The reconstruction of a τ h lepton starts with a PF jet as the initial seed. This is followed by the reconstruction of the π 0 components in the jet which are then combined with the charged hadron components to fully reconstruct the decay mode of the τ h and to calculate its fourmomentum [59]. The identification of τ h is performed by a MVA based discriminant [60]. The main handle to separate hadronic τ decays from quark and gluon jets is the isolation of the τ h candidate from other charged hadrons and photons. Variables that are sensitive to the distance of separation between the production and decay vertices of the τ h candidate complement the MVA inputs. This algorithm achieves a τ h identification efficiency of 50% with a misidentification rate for quark and gluon jets below 1%. Additional discriminants are used to separate τ h candidates from electrons and muons [60]. The discriminant against electrons uses variables sensitive to electron shower shape, electron track, and τ h decay kinematics. The discriminant against muons uses inputs based on calorimetric information of the τ h jet and reconstructed hits and track segments in the muon system.

HH mass reconstruction and event selection
This analysis is based on data satisfying a τ h τ h trigger which requires the presence of two τ h objects with a p T threshold of 35 GeV and η ≤ 2.1 for each τ h . A further selection of events is made offline. It is first ensured that the data considered in the analysis are of good quality and each event contains a primary vertex with the absolute value of the z coordinate less than 24 cm, and within the radial distance of 2 cm from the beam axis. The following analysis specific selection criteria are then applied, determined by the need to suppress specific types of backgrounds. These selection criteria depend on the mass of the pair of τ h candidates and the pair of b quark jets which are determined as follows.
The H that decays into a pair of τ h leptons is reconstructed by a likelihood based algorithm, referred to as SVfit [61]. The algorithm uses the four-momenta of the two τ h candidates, the magnitude and direction of the p miss T vector as well as the event-by-event estimate of the p miss T resolution as input to reconstruct the full four-momentum vector (p T , η, φ, and mass) of the pair of τ h candidates without any constraint on its mass. A mass window constraint is later applied as described below. The four-vector of the H that decays into b quarks is reconstructed by means of a kinematic fit. The fit varies the energy of the highest quality (according to the CSV algorithm) b quark jet within the expected resolution, keeping the jet direction fixed, subject to the constraint that the invariant mass of the two b quark jets equals m H = 125 GeV. Further selection is based on a mass window criterion as described below.
In the search for resonant HH production, the four-momentum vectors of the two H's are used to reconstruct the mass of the HH system, m HH . We assume that the width of the new particle X is small compared to the experimental resolution on the mass of the H pair, which, for resonances of true mass m X in the range 300 GeV to 1000 GeV, typically amounts to 8% times m X . A peak in the HH mass distribution is expected this case. The search for heavy spin 0 and spin 2 resonances is hence based on finding a peak in the HH mass spectrum.
In the non-resonant case, the mass distribution of the H pair is expected to be broader than the experimental resolution. After comparing different observables in terms of their capability to separate a potential signal from the background we have found that the observable m T2 [62] performs the best. Our search for non-resonant HH production is hence based on the m T2 variable which is an analog of the transverse mass variable used in W → ν analyses, adapted to the cascade decays of tt pairs to pairs of b quarks, leptons, and neutrinos. It improves the separation of the HH signal in particular from the tt background, due to the fact that values of the m T2 variable extend up to 300-400 GeV for signal events, while for tt background events they are concentrated below the top quark mass. The usage of this observable in analyses of non-resonant HH production in the bbττ final state was first proposed in Ref. [63].

Definition of event categories
The selection of events is based on the following additional requirements: • The event is required to contain two τ h candidates with p T > 45 GeV and |η| < 2.1, which pass the identification criteria described in Section 3.2. Both τ h candidates are required to be matched to the τ objects that trigger the event within ∆R < 0.5.
Here ∆R = √ (∆η) 2 + (∆φ) 2 and ∆η and ∆φ are the distances in pseudorapidity and azimuthal angle (in radians), respectively, between the reconstructed tau object and the tau object at the trigger level.
• The two τ h candidates are required to be of opposite charge. The τ h τ h invariant mass (m ττ ), reconstructed by the SVfit algorithm, is required to be in the window 80-140 GeV. If multiple combinations exist in an event, the combination with the highest sum of outputs from the MVA based discriminant that separates the τ h candidate from quark and gluon jets, is taken.
• The event is required to contain two jets of p T > 20 GeV and |η| < 2.4. The jets are required to be separated from each of the two τ h candidates by ∆R > 0.5. The mass of the two jets is required to be within the window 80 < m jj < 170 GeV. • Events containing an isolated electron of p T > 15 GeV and |η| < 2.4, or an isolated muon of p T > 15 GeV and |η| < 2.4 are rejected.
In the search for non-resonant HH production, the Lorentz boost of the H's and the resulting boost of the τ h lepton pair coming from their decays is used to further distinguish between signal and background events by requiring the distance in η-φ between the two τ h candidates, ∆R ττ , to be less than 2.0. This criterion is not used in the resonant HH search in order to preserve sensitivity in the low mass (m HH < 500 GeV) region. Except for the ∆R ττ criterion, the event and object selection applied in the search for non-resonant and for resonant HH production are identical.

Definition of event categories
The HH → bbττ signal events are expected to contain two b quark jets in the final state. The efficiency to reconstruct a single b jet is higher than reconstructing two b jets in an event. The efficiency of signal selection is therefore enhanced in this analysis by accepting events with one b tagged jet and one jet which is not b tagged. A control region containing events with two or more jets, none of which passes the b tagging criteria, is used to constrain systematic uncertainties. More specifically, the event categories are: Events in this category are required to contain at least two jets of p T > 20 GeV and |η| < 2.4 which are selected by the CSV discriminant described in Section 3.1.
Events in this category are required to contain one jet of p T > 20 GeV and |η| < 2.4, which is selected by the CSV discriminant and one or more additional jets of p T > 20 GeV. These jets are required to either not satisfy |η| < 2.4 or not to be selected by the CSV discriminant.
• 0 b tags Events in this category are required to contain at least two jets of p T > 20 GeV, all of which either do not satisfy |η| < 2.4 or are not selected by the CSV discriminant.
These categories are mutually exclusive. For the purpose of studying the modeling of data by MC simulation in a region that is not sensitive to the presence or the absence of signal events, we define as 'inclusive' category the union of all three categories. No selection criteria are applied on m ττ , m jj , or ∆R ττ in the inclusive category.

Background estimation
The two important sources of background in the 0 b tag and 1 b tag categories are events containing Z/γ * → ττ decays and multijet production. In the 2 b tag category Z/γ * → ττ decays and tt events are dominant sources of background events.

The multijet events
The reconstructed τ h candidates in multijet events are typically due to the misidentification of quark or gluon jets. The contribution from this background in the signal region, in terms of event yield and shape of the distributions in m HH and m T2 ("shape template"), is determined entirely from data. The normalization and shape is obtained separately in each event category, from events that pass the selection criteria described in Sections 4 and contain two τ h candidates of opposite charge. It is required that the leading (higher p T ) τ h candidate passes relaxed, but fails the nominal τ h identification criteria. The probabilities for the leading τ h candidate to pass the relaxed and nominal τ h identification criteria are measured in events that contain two τ h candidates of the same charge, as functions of p T of the leading τ h candidate in three regions of η, |η| < 1.2, 1.2 < |η| < 1.7, and 1.7 < |η| < 2.1. A linear function is fitted to the variation of the ratio of these two probabilities with p T and is applied as an event weight to obtain the estimate for the shape template of the multijet background in the signal region. Contributions from other backgrounds to these events are subtracted based on MC predictions.

The Z/γ * → ττ events
The dominant irreducible Z/γ * → ττ background in the event categories with 2 b tags, 1 b tag, and 0 b tags is modeled by applying embedding to Z/γ * → µµ events selected from data as described in Section 2.2. The embedded sample is normalized to the Z/γ * → ττ event yield obtained from the MC simulation in the inclusive event category. The correction due to tt contamination is performed by subtracting the distribution in m HH or m T2 whose shape and normalization are determined using the tt embedded sample from that in the Z/γ * → ττ embedded sample in each event category. An uncertainty on the number of events in each bin is set to the sum of uncertainties of the Z/γ * → ττ and tt embedded yields in that bin, added in quadrature.
The embedded samples cover only a part of the Z/γ * → ττ background, namely events in which both reconstructed τ h candidates match generator-level hadronic τ decays, because of requirements that are applied at the generator level during the production of the embedded samples to enhance the number of events that pass the selection criteria described in Sections 4. The small additional contribution arising from Z/γ * → ττ production in which one or both reconstructed τ h candidates are due to a misidentified electron, muon, or jet are taken from the Z/γ * → ττ MC sample.

Other backgrounds
The contribution of tt background is estimated using an MC sample after reweighting the events as described in Section 2.2. The background contributions arising from W+jets, Z/γ * → ( = e, µ), single top quark, and di-boson production, as well as from the production of events with a single SM H boson are small and are modeled using MC samples.

Systematic uncertainties Systematic uncertainties
The systematic uncertainties in this analysis may affect the number of signal or background events selected in a given event category or affect the relative number of signal or background events in individual bins of kinematic distributions. An additional uncertainty arises due to the limited statistics available to model the m HH or m T2 distributions of individual backgrounds in some of the event categories. The treatment of such uncertainties is described in Section 8. The systematic uncertainties relevant to this analysis are • The τ h trigger and identification efficiency The uncertainty in the τ h identification efficiency has been measured as 6% using Z/γ * → ττ → µτ h events. The τ h candidates in Z/γ * → ττ events typically have p T in the range 20 to 50 GeV. An uncorrelated uncertainty of 20%p T /(1000 GeV) is added to account for the extrapolation to the high-p T region, including the uncertainty in the charge misidentification rate of high-p T τ leptons. The above uncertainties have been taken from Ref. [60]. The uncertainty in the efficiency of the τ h τ h trigger amounts to 4.5% per τ h candidate [24].
• τ h energy scale The uncertainty in the τ h energy scale is taken as 3% [60].
• Background yields The rate of the Z/γ * → ( = e, µ) background is attributed an uncertainty of 5%. The normalization of the Z/γ * → ττ embedded samples, as described in Section 6.2, is attributed an uncertainty of 5%. An additional uncertainty of 5% is assigned to the fraction of Z/γ * → ττ events entering the 2 b tags and 1 b tag categories. This uncertainty has been introduced to cover potential small biases of the embedding technique. The rate of the tt background is known with an uncertainty of 7%. The uncertainty in the MC yield of single top quark and di-boson backgrounds amounts to 15%. An uncertainty of 30% has been applied to the W+jets background yield obtained from MC. The above uncertainties have been taken from Refs. [24, 64].
• Integrated luminosity The uncertainty in the integrated luminosity is taken as 2.6% [65]. This uncertainty is applied to signal and to Z/γ * → ( = e, µ, τ), W+jets, single top quark and di-boson backgrounds. This uncertainty is not applied to the tt background, as this background is normalized to the top quark pair production cross section measured by CMS with a correction factor obtained from a tt dominated control region in data as described in Section 2.2. The normalization of the multijet background is obtained from data and hence is not subject to the luminosity uncertainty.

• Jet energy scale
Jet energy scale uncertainties range from 1 to 10% and are parametrized as functions of jet p T and η [55]. They affect the yield of signal and background events in different event categories and the shape of the m HH and m T2 distributions.
• The b tagging efficiency and the mistag rate Uncertainties in the b tagging efficiencies and the mistag rates result in event migration between categories. These are evaluated as functions of jet p T and η as determined in Ref. [56] and are applied to MC samples.
• The multijet background estimation The uncertainty in this background contribution is obtained by adding the statistical uncertainty in the yield of events in the sample with two opposite charge τ h candidates in quadrature with the uncertainty in the slope and offset parameters of the function used as event weight to the shape template as described in Section 6.1.
• The p miss T resolution and response The uncertainties related to the magnitude and direction of the p miss T vector, which affect the shape of the m HH and m T2 distributions, are covered by uncertainties in the Z boson recoil correction. The Z boson recoil correction is computed by comparing data with simulation in Z → ee, Z → µµ, and photon+jets samples, which do not have any genuine missing transverse momentum. All observables related to p miss T (including m HH and m T2 ) are recomputed by varying p miss T within its uncertainty [33] and applied to MC samples.

• The top quark p T reweighting
The reweighting that is applied to simulated tt events (Section 2.2) is varied between one (no correction) and twice the reweighting factor (overcorrection by 100%) to account for the uncertainty due to reweighting [39,40].
• Other sources The uncertainties on the SM HH cross section are +4.1%/−5.7% due to scale, ±5% due to approximations concerning top quark mass effects that are made in the theoretical calculations, ±2.6% due to α S and ±3.1% due to the parton density function [3]. The uncertainty due to the H → ττ (H → bb) branching fraction is ±3.3% (±3.2%) [66]. The effect of the uncertainty on the number of pileup interactions amounts to less than 1% and is neglected.

Signal extraction
Signal rates are determined from a binned maximum likelihood fit for signal plus background and background-only hypotheses. In case of resonant (non-resonant) HH production, we fit the distribution of m HH (m T2 ), reconstructed as described in Section 4. Constraints on systematic uncertainties that correspond to multiplicative factors on the signal or the background yield (e.g. cross sections, efficiencies, misreconstruction rates, and sideband extrapolation factors) are represented by log-normal probability density functions. Systematic uncertainties in the shape of m HH and m T2 distributions for signal as well as background processes are accounted for by the 'vertical template morphing' technique [67] and represented by Gaussian probability density functions. The Barlow-Beeston method [67,68] is employed to account for statistical uncertainties on the m HH and m T2 shape templates.

Observed yields
The number of events observed in the event categories with 2 b tags, 1 b tag, and 0 b tags as well as the expected yield of background processes in these categories are given in Table 1. The signal rate expected for non-resonant HH production has been computed for a cross section σ(pp → HH) of 1 pb, corresponding to 100 times the SM cross section, and SM event kinematics [69,70]. In the case of resonant HH production, the signal yield has been computed for a resonance X (radion or graviton) of mass m X = 500 GeV and a σ(pp → X For non-resonant HH production the distributions of m T2 are shown in Fig. 3. For the resonant case the distribution of m HH for events selected in the three categories mentioned above are shown in Fig. 4. In both figures, the sum of W+jets, single top quark and di-boson events and of Z+jets events in which one or both reconstructed τ h are due to a misidentified e, µ, or jet is referred to as "electroweak" background. Bins in which zero events are observed in the data are indicated by the absence of a data point. The vertical bar drawn in these bins indicate the 84% confidence interval, corresponding to a tail probability of 16%. The event yields and the shape of mass distributions observed in data are in agreement with background predictions.
No evidence for the presence of a signal is observed.

Cross section limits
We have set 95% CL upper limits on cross section times branching fraction for HH production using a modified frequentist approach, known as the CL s method [71][72][73]. For non-resonant production SM event kinematics have been assumed. Some model dependency is expected in this case, as the signal acceptance times efficiency as well as the shape of the m T2 distribution vary as functions of the m HH spectrum predicted by the model. The observed (expected) limits on σ(pp → HH) are 0.59 (0.94 +0.46 −0.24 ) pb, corresponding to a factor of about 59 (94) times the cross section predicted by the SM. For the production of resonances decaying to a pair of SMlike H's of mass m H = 125 GeV the difference between the limits computed for radion → HH and graviton → HH signals is small, indicating that the limits on resonant HH production cross section do not depend on these particular models. The limits obtained for resonant HH production are given in Table 2 and are shown in Fig. 5. In this figure, the expected limits are computed for a generic spin 0/2 resonance decaying to two SM H's. The theoretical curves for the graviton case are based on KK graviton production in the bulk and RS1 models, respectively [18,19]. To obtain the radion theoretical curves, cross section for radion production via gluon fusion are computed (to NLO electroweak and NNLO QCD accuracy) for different values of the fundamental theoretical parameter Λ R . These values are then multiplied by a k factor calculated for SM-like H production through gluon-gluon fusion [74][75][76]. Table 2: The 95% CL upper limits on resonant HH production (σ(pp → X) B(X → HH)) in units of pb for spin 0 (radion) and spin 2 (graviton) resonances X, at different masses m X , obtained from the HH search in the decay channel bbττ. The results of the search for HH production in the bbττ decay channel are combined with those in the decay channels γγbb and bbbb, published in Refs. [25,26] respectively. The combination is performed by adding the three individual log likelihood functions. The correlated systematics are taken into account by using the same nuisance parameters for the fully correlated sources. They are the luminosity uncertainty, the uncertainty on the b tagging efficiency, the uncertainties related to the underlying event and parton showering, the uncertainties on the   branching fractions of the three HH decays channels, and the theoretical uncertainties on the SM non-resonant HH cross section, parton density functions and α S . The uncertainty on the branching fraction of H → γγ is ±5% [66].
The signal yield in the three decay channels is determined assuming that the branching fractions for the decays H → bb, H → ττ, and H → γγ are equal to the SM predictions [66] for a H with mass m H = 125 GeV. The datasets analyzed by the γγbb and bbbb decay channels correspond to integrated luminosities of 19.7 and 17.9 fb −1 , recorded at √ s = 8 TeV respectively. The search in the γγbb decay channel targets resonant as well as non-resonant HH production, while the search in the bbbb decay channel focuses on resonant HH signals. No evidence for a signal is observed in the combined search.
The limits on resonant HH production obtained from the combination of bbττ, γγbb, and bbbb channels are given in Table 3 and Fig. 6. In the case of non-resonant HH production, an observed (expected) limit on σ(pp → HH) of 0.43 pb (0.47 +0.20 −0.12 pb), corresponding to 43 (47) times the SM cross section, is obtained by combining the bbττ and γγbb decay channels. The low mass sensitivity (m HH ≤ 400 GeV) is dominated by the γγbb channel while the high mass (m HH > 700 GeV) sensitivity is driven by the bbbb channel. The bbττ channel is competitive with the γγbb channel in the intermediate mass range (400 GeV < m HH ≤ 700 GeV).

Summary
A search has been performed for events containing a pair of SM-like H's in resonant and nonresonant production of the pair in the channel where one boson decays to a pair of b quarks and the other to a τ lepton pair, in pp collisions collected by the CMS experiment at 8 TeV center-of-mass energy, corresponding to an integrated luminosity of 18.3 fb −1 . Results are ex- Table 3: The 95% CL upper limits on resonant HH production (σ(pp → X) B(X → HH)) in units of fb for spin 0 (radion) and spin 2 (graviton) resonances X, at different masses m X , obtained from the combination of HH searches performed in the bbττ, γγbb, and bbbb decay channels.   Figure 6: 95% CL observed and expected upper limits on the cross section times branching fraction (σ(pp → X) B(X → HH)) for a spin 0 (left) and for a spin 2 (right) resonance X as functions of the resonance mass m X , obtained from the combination of searches performed in the bbττ, γγbb and bbbb decay channels. The green and yellow bands represent, respectively, the 1 and 2 standard deviation extensions beyond the expected limit. Also shown are theoretical predictions corresponding to WED models for radions for values of Λ R = 1, 3 and for RS1 and Bulk KK gravitons [18,19]. The other WED model parameters are kl = 35 and k/M Pl = 0.2, assuming an elementary top hypothesis and no radion-Higgs (r/H) mixing. pressed as 95% CL upper limits on the production of a signal. The limit on non-resonant HH production corresponds to a factor of 59 times the rate expected in the SM. For resonant X → HH production, the limit on σ(pp → X) B(X → HH) for a resonance of spin 0 and spin 2 ranges, respectively, from 5.42 and 3.97 pb at a mass m X = 300 GeV to 0.14 pb and 0.14 pb at m X = 1000 GeV.
The results of the search in the bbττ decay channel are combined with those in the γγbb and bbbb decay channels. For non-resonant HH production, the combination of bbττ and γγbb decay channels yields a limit that is a factor of 43 times the SM rate. The limit on resonant HH production obtained from the combination ranges from 1.13 and 1.09 pb at m X = 300 GeV, to 21 and 18 fb at m X = 1000 GeV for resonances of spin 0 and spin 2 respectively.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Aus  [23] CMS Collaboration, "Searches for heavy Higgs bosons in two-Higgs-doublet models and for t→ ch decay using multilepton and diphoton final states in pp collisions at 8