Search for Higgs boson pair production in events with two bottom quarks and two tau leptons in proton-proton collisions at sqrt(s) = 13 TeV

A search for the production of Higgs boson pairs in proton-proton collisions at a centre-of-mass energy of 13 TeV is presented, using a data sample corresponding to an integrated luminosity of 35.9 inverse femtobarns collected with the CMS detector at the LHC. Events with one Higgs boson decaying into two bottom quarks and the other decaying into two tau leptons are explored to investigate both resonant and nonresonant production mechanisms. The data are found to be consistent, within uncertainties, with the standard model background predictions. For resonant production, upper limits at the 95% confidence level are set on the production cross section for Higgs boson pairs as a function of the hypothesized resonance mass and are interpreted in the context of the minimal supersymmetric standard model. For nonresonant production, upper limits on the production cross section constrain the parameter space for anomalous Higgs boson couplings. The observed (expected) upper limit at 95% confidence level corresponds to about 30 (25) times the prediction of the standard model.


Introduction
The discovery of the Higgs boson (H) by the ATLAS and CMS Collaborations [1][2][3] was a major step towards improving the understanding of the mechanism of electroweak symmetry breaking (EWSB).With the mass of the Higgs boson now precisely determined [4], the structure of the Higgs scalar field potential and the Higgs boson self-couplings are precisely predicted in the standard model (SM).While the measured properties of the Higgs boson are thus far consistent with the expectations from the SM [5], the measurement of the Higgs boson self-coupling provides an independent test of the SM and verification that the Higgs mechanism is truly responsible for the EWSB by giving access to the shape of the Higgs scalar field potential [6].
The trilinear self-coupling of the Higgs boson (λ HHH ) can be extracted from the measurement of the Higgs boson pair (HH) production cross section.In the SM, for proton-proton (pp) collisions at the CERN LHC, this process occurs mainly via gluon-gluon fusion and involves either couplings of the Higgs boson to virtual fermions in a quantum loop, or the λ HHH coupling itself, with the two processes interfering destructively as illustrated in Fig. 1.
The SM prediction for the cross section is σ HH = 33.49+4.3% −6.0% (scale) ± 5.9% (theo) fb [7][8][9][10][11].This value was computed at the next-to-next-to-leading order (NNLO) of the theoretical perturbative quantum chromodynamics (QCD) calculation, including next-to-next-to-leadinglogarithm (NNLL) corrections and finite top quark mass effects at next-to-leading order (NLO).The theoretical uncertainties in σ HH include uncertainties in the QCD factorization and renormalization scales, the strong coupling parameter α S , parton distribution functions (PDF), and unknown effects from the finite top quark mass at NNLO. 1: Cartoon of the region in the plane (g ⇤ , /g ⇤ ), defined by Eqs. ( 13), (14), that can be probed y an analysis including only dimension-6 operators (in white).No sensible e↵ective field theory escription is possible in the gray area ( < g min ), while exploration of the light blue region g min < < p g ⇤ g min ) requires including the dimension-8 operators.erivative terms (which correspond to dimension-8 operators in the limit of linearly-realized W symmetry).The e↵ect of the neglected derivative operators will be then studied by nalyzing their impact on angular di↵erential distributions and shown to be small in our ase due to the limited sensitivity on the high m hh region.
The Feynman diagrams that contribute to the gg !hh process are shown in Fig. 1: Cartoon of the region in the plane (g ⇤ , /g ⇤ ), defined by by an analysis including only dimension-6 operators (in white) description is possible in the gray area ( < g min ), while ex (g min < < p g ⇤ g min ) requires including the dimension-8 opera derivative terms (which correspond to dimension-8 operator EW symmetry).The e↵ect of the neglected derivative op analyzing their impact on angular di↵erential distribution case due to the limited sensitivity on the high m hh region.
The Feynman diagrams that contribute to the gg !hh p diagram is characterized by a di↵erent scaling at large ene 1: Cartoon of the region in the plane (g ⇤ , /g ⇤ ), defined by Eqs. ( 13), (14), that can be probed y an analysis including only dimension-6 operators (in white).No sensible e↵ective field theory escription is possible in the gray area ( < g min ), while exploration of the light blue region g min < < p g ⇤ g min ) requires including the dimension-8 operators.erivative terms (which correspond to dimension-8 operators in the limit of linearly-realized W symmetry).The e↵ect of the neglected derivative operators will be then studied by nalyzing their impact on angular di↵erential distributions and shown to be small in our ase due to the limited sensitivity on the high m hh region.
The Feynman diagrams that contribute to the gg !hh process are shown in Fig.Beyond the standard model (BSM) physics effects can appear either via anomalous couplings of the Higgs boson or via new particles that can be directly produced or contribute to the quantum loops responsible for HH production.The experimental signature would be an enhancement of the HH production cross section for a specific value of the invariant mass of the pair (resonant production) or over the whole invariant mass spectrum (nonresonant production).
Resonant double Higgs boson production is predicted by many extensions of the SM such as the singlet model [12][13][14], the two-Higgs-doublet model [15] and its realisation as the minimal supersymmetric standard model (MSSM) [16,17], and models with warped extra dimensions (WED) [18,19].Although the physics motivation and the phenomenology of these theoretical models are very different, the signal is represented by a CP-even scalar particle (S) decaying into a Higgs boson pair, with an intrinsic width that is often negligible with respect to the detector resolution.
In the nonresonant case, the BSM physics is modelled through an effective Lagrangian that extends the SM Lagrangian with dimension-6 operators [20].Five Higgs boson couplings result from this parametrization: the Higgs boson coupling to the top quark, y t , the trilinear coupling λ HHH , and three additional couplings, denoted as c 2 , c 2g , and c g using the notation in Ref. [7], that represent, respectively, the interactions of a top quark pair with a Higgs boson pair, of a gluon pair with a Higgs boson pair, and of a gluon pair with a single Higgs boson.For simplicity, we investigate only anomalous y t and λ HHH couplings, while the other anomalous couplings are assumed to be zero, and parametrize the deviations from the SM values as k λ = λ HHH /λ SM HHH and k t = y t /y SM t .Extension of these results to any combination of the couplings can be obtained by following the procedure detailed in Ref. [21].These two couplings are currently largely unconstrained by experimental results, and deviations from the SM can be accommodated by the combined measurements of Higgs boson properties [5] depending on the particular assumptions made about the BSM physics contributions.
Previous searches for the production of Higgs boson pairs were performed by both the AT-LAS [22,23] and CMS [24,25] Collaborations using the LHC data collected at √ s = 8 and 13 TeV.The most sensitive upper limit at 95% confidence level (CL) on HH production corresponds to 43 times the rate predicted by the SM and is obtained from the combination of the HH → bbγγ and HH → bbτ + τ − decay channels using data collected at √ s = 8 TeV [26].
In this Letter we present a search for Higgs boson pair production in the final state where one Higgs boson decays to bb and the other decays to τ + τ − .For simplicity, we refer to this process as HH → bbττ in the following, omitting the quark and lepton charges.This process has a combined branching fraction of 7.3% for a Higgs boson mass of 125 GeV.Its sizeable branching fraction, together with the relatively small background contribution from other SM processes, makes this final state one of the most sensitive to HH production.Three final states of the τ lepton pair are considered: one of the two τ leptons is required to decay into hadrons and a neutrino (τ h ), while the other can decay either to the same final state, or into an electron (τ e ) or a muon (τ µ ) and neutrinos.Together, these three final states include about 88% of the decays of the ττ system and are the most sensitive ones for this search.The data sample analyzed corresponds to an integrated luminosity of 35.9 fb −1 collected in pp collisions at √ s = 13 TeV.
The search described in this Letter improves on the previous HH → bbττ results [26] by including final states with a leptonic τ decay, improving the event categorization, introducing multivariate methods for the background rejection, and optimizing the event and object selection for the LHC collisions at √ s = 13 TeV.

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 solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections.Forward calorimeters extend the pseudorapidity coverage provided by the barrel and endcap detectors.Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.Events of interest are selected using a two-tiered trigger system [27].The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a time interval of less than 4 µs.The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to less than 1 kHz before data storage.A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, including pseudorapidity η and azimuthal angle ϕ, can be found in Ref. [28].

Modelling of physics processes
Simulated samples of resonant and nonresonant HH production via gluon-gluon fusion are generated at leading order (LO) precision with MADGRAPH5 aMC@NLO 2.3.2 [29].In the case of resonant production, separate samples are generated for mass values of the resonance ranging from 250 to 900 GeV.In the case of nonresonant production, separate samples are generated for different values of the effective Lagrangian couplings, including the couplings predicted by the SM [21,30].In the latter case, an event weight determined as a function of the generated HH pair kinematics is applied to these samples to model signals corresponding to additional points in the effective Lagrangian parametrization.
Backgrounds arising from Z/γ * → and W → ν in association with jets (with = e, µ, τ), diboson (WW, ZZ, and WZ), and SM single Higgs boson production are simulated with MAD-GRAPH5 aMC@NLO 2.3.2 at LO with MLM merging [31], while the single top and tt backgrounds are simulated at NLO precision with POWHEG 2.0 [32,33].The NNPDF3.0 [34] PDF set is used.In order to increase the number of simulated events that satisfy the requirements detailed in Section 4, the inclusive simulation of the Z/γ * and W processes is complemented by samples simulated in selected regions of multiplicity, flavour, and the transverse momentum scalar sum of the partons emitted at the matrix element level.Signal and background generators are interfaced with PYTHIA 8.212 [35] with the tune CUETP8M1 [36] to simulate the multiparton, parton shower, and hadronization effects.The simulated events include multiple overlapping hadron interactions as observed in the data.
The tt, Z/γ * → , W → ν and single top quark samples are normalized to their theoretical cross sections at NNLO precision [37][38][39], and the diboson samples are normalized to their cross section at NLO precision [40].The single Higgs boson production cross section is computed at the NNLO precision of the QCD corrections and at the NLO precision of electroweak corrections [7,[41][42][43][44].

Object reconstruction and event selection
In order to reconstruct an HH → bbττ candidate event, it is necessary to identify the e, µ, and τ h leptons, the jets originating from the two b quarks, and the missing transverse momentum vector p miss T , defined as the projection onto the plane perpendicular to the beam axis of the negative vector sum of the momenta of all reconstructed particle-flow objects in an event.Its magnitude is referred to as p miss T .The particle-flow (PF) event algorithm [45] reconstructs and identifies each individual particle (PF candidate) with an optimized combination of information from the various elements of the CMS detector.The momentum of the muons is obtained from the curvature of the corresponding track.The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex, as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track.The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers.Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.Complex objects, such as τ h , jets, and the p miss T vector, are reconstructed from PF candidates.For each event, hadronic jets are clustered from PF candidates with the infrared and collinear safe anti-k T algorithm [46,47], operated with distance parameters of 0.4 and 0.8.These jets are denoted as "AK4" and "AK8" in the following.Leptons from b hadron decays within a jet are considered as constituents by the algorithm.The jet momentum is determined as the vectorial sum of all particle momenta in the jet, and is found in the simulation to be within 5 to 10% of the true momentum over the whole p T spectrum and detector acceptance.The invariant mass of AK8 jets is obtained by applying the soft drop jet grooming algorithm [48,49], that iteratively decomposes the jet into subjets to remove the soft wide-angle radiation and mitigates the contribution from initial state radiation, underlying event, and multiple hadron scattering.Jet energy corrections are derived from the simulation, and are confirmed with in situ measurements using the energy balance of dijet, multijet, γ+jet, and leptonic Z+jet events [50,51].The PF components of the jets are used to reconstruct τ h candidates using the hadrons plus strips algorithm [52,53], combining either one or three charged particle tracks with clusters of photons and electrons to identify the decay mode of the τ lepton.
Events in the bbτ µ τ h (bbτ e τ h ) final state have been recorded using a set of triggers that require the presence of a single muon (electron) in the event.The selected events are required to contain a reconstructed muon (electron) [54,55] of p T > 23 (27) GeV and |η| < 2.1 and a reconstructed τ h candidate [52] of p T > 20 GeV and |η| < 2.3.The muon (electron) candidate must satisfy the relative isolation requirement I rel < 0.15 (0.1) [54,55], while the τ h candidate must satisfy the "medium" working point of a multivariate isolation discriminant [52], that corresponds to a signal efficiency of about 60% and a jet misidentification rate ranging between 0.1% and 1% depending on the jet p T .The reconstructed tracks associated to the selected electron, muon, and τ h candidates must be compatible with the primary pp interaction vertex of the event.Electrons and muons erroneously reconstructed as a τ h candidate are rejected using discriminants based on the information from the calorimeters and muon detectors and on the properties of the PF candidates that form the τ h candidate, as is detailed in [52].
A trigger requiring the presence of two τ h candidates is used to record events in the bbτ h τ h final state.The selected events must contain two reconstructed τ h candidates with p T > 45 GeV and |η| < 2.1, that are required to pass the "medium" working point of the multivariate isolation discriminant and whose associated tracks must be compatible with the primary pp interaction vertex of the event.The discriminants that suppress the contribution from prompt electrons and muons are applied to both τ h candidates as in the bbτ µ τ h and bbτ e τ h final states.
For all three final states, the two selected τ leptons are required to have opposite electric charge.Events containing additional isolated muons or electrons are rejected to reduce the Z/γ * → background contribution.
Events selected with the criteria described above (τ µ τ h , τ e τ h , τ h τ h ) are required to have two additional AK4 jets with p T > 20 GeV and |η| < 2.4.In the case of HH production via a resonance of mass 700 GeV or higher, the two jets originating from the H → bb decay partially overlap due to the high Lorentz boost of the Higgs boson, and are reconstructed at the same time as two separate AK4 jets and as a single AK8 jet.To profit from this information, the event is classified as "boosted" if it contains at least one AK8 jet of invariant mass larger than 30 GeV and p T > 170 GeV that is composed of two subjets, each geometrically matched to one of the selected AK4 jets (∆R(AK4, subjet) < 0.4, where ∆R = √ (∆η) 2 + (∆ϕ) 2 denotes the spatial separation of the jet candidates).The event is classified as "resolved" if any of these requirements is not satisfied.This classification provides a clear separation of the signal topology against the tt background, where the two jets are typically more spatially separated and not reconstructed as a single AK8 jet.The AK8 jet mass requirement is applied to reject candidates resulting from a single quark or gluon hadronization or poorly reconstructed by the soft drop algorithm.
The combined secondary vertex [56] algorithm is applied to the selected jets to identify those originating from a bottom quark and reduce the contribution from the multijet background where jets are initiated by light quarks or gluon radiation.Both the "medium" and the "loose" working points of the b tagging discriminant [57] are used in this search as described below.The efficiency and rate of erroneous b jet identification are about 60% (80%) and 1% (10%) respectively for the "medium" ("loose") working point.
Jets reconstructed in events classified as "resolved" are defined as b-tagged if they satisfy the "medium" working point of the b tagging algorithm.These events are classified into two groups according to the number of b-tagged jets: the group with at least two b-tagged jets (2b) has the best sensitivity, and the group with exactly one b-tagged jet (1b1j) increases the signal acceptance.Both AK4 jets previously selected in the events classified as "boosted" are required to satisfy the "loose" working point of the b tagging discriminant.

Signal regions and discriminating observables
After the object selection and event classification, the kinematic information of the event is exploited to reduce the contribution from background processes.The invariant mass of the two τ lepton candidates, m ττ , is reconstructed using a dynamic likelihood technique called SVfit [58] that combines the kinematics of the two visible lepton candidates and the missing transverse momentum in the event.The bb invariant mass, m bb , is estimated from the two selected jet candidates for "resolved" topologies and from the invariant mass of the AK8 jet for "boosted" topologies.In the "resolved" case, the events are required to satisfy the condition: where the values of 35 and 45 GeV are related to the mass resolution of the ττ and bb systems and 116 and 111 GeV correspond to the position of the expected reconstructed 125 GeV Higgs boson peak in the m ττ and m bb distributions, respectively.The selection has been optimized for the SM HH process to obtain a signal efficiency of approximately 80% and a background reduction of about 85% in the most sensitive event categories.The m bb peak is shifted below the Higgs boson mass value because the momenta of neutrinos from b hadron decays are not measured.This effect also prevents the SVfit algorithm from fully recovering the ττ system mass value.In the "boosted" case the events are required to satisfy: In addition to the previous requirements, a multivariate discriminant is applied to the events in the resolved categories of the bbτ µ τ h and bbτ e τ h final states to identify and reject the tt process, which is the most important source of background.The discriminant is built using the boosted decision tree (BDT) [59,60] algorithm that is trained on a combination of τ µ τ h and τ e τ h simulated signal and background events.The algorithm identifies the kinematic differences between the two processes and assigns to every selected event a number that defines its compatibility with a signal or background topology.Two separate BDT trainings are performed to achieve an optimal performance for all the signal processes studied.
One training is performed using resonant signals with masses m S ≤ 350 GeV as input.Eight variables are used in the discriminant training because of their good separation between ) 2 denotes the transverse mass of the selected lepton candidate, with a similar definition for m T (τ h ).The ∆R separations of the two b quarks and of the two tau leptons are multiplied by the H bb and H ττ candidate p T respectively to reduce their dependence on the m S hypothesis.All the selected variables contribute significantly to the discrimination achieved with the trained BDT.The same training is used both for the search for resonant HH production up to m S = 350 GeV and for the search for nonresonant HH production.No loss of performance is observed by using this training in comparison to a dedicated training on nonresonant signals.Different selections on the BDT discriminant output are applied in the two searches to maximize the sensitivity: these selections correspond to a rejection of the tt background of approximately 90 and 70% for the resonant and nonresonant searches, respectively, for a signal efficiency ranging between 65 and 95% depending on the signal hypothesis considered.
A second training is performed on the resonant signals of mass m S > 350 GeV.The variables used as inputs to this training are the same as in the previous case, but replacing ∆R(b, b) p T (H bb ) and ∆R( , τ h ) p T (H ττ ) with ∆R(b, b) and ∆R( , τ h ).The selection on the BDT output is chosen to maximize the sensitivity and corresponds to a rejection of the tt background of approximately 90% for a signal efficiency ranging between 70 and 95% depending on the value of m S .In the case of the resonant search, the selections applied to the two BDT discriminants define low-mass (LM) and high-mass (HM) signal regions.
In the resonant search, the invariant mass of the two visible τ lepton decay products and the two selected b jets is used to search for a possible signal above the expected background event distribution.In order to improve the resolution and to enhance the sensitivity of the analysis, the invariant mass is reconstructed using a kinematic fit (m KinFit HH ) that is detailed in Ref. [61].The fit is based on the four-momenta of the τ and b candidates and on the p miss T vector in the event, and is performed under the hypothesis of two 125 GeV Higgs bosons decaying into a bottom quark pair and a τ lepton pair.The use of the kinematic fit improves the resolution on m HH by about a factor of two compared to the four-body invariant mass of the reconstructed leptons and jets.
The stransverse mass or m T2 variable is used in the search for a nonresonant signal.This variable, originally introduced for supersymmetry searches involving invisible particles in the final state [62,63] and later proposed for HH searches in bbττ events [64], is used to reconstruct events where two equal mass particles are produced and each undergoes a two-body decay into a visible and an invisible particle.The m T2 variable is defined as the largest mass of the parent particle that is compatible with the kinematic constraints of the event.In the case of the bbττ decay, where the dominant background is tt production, the parent particle is interpreted as the top quark that decays into a bottom quark and a W boson.Following the description in Ref. [64], we denote with b, b the momenta of the two selected b jets and with m b , m b their invariant masses, and we introduce the c, c symbols to denote the momenta of the other particles produced in the top quark decay corresponding to the measured leptons and the neutrinos.We also set m c = m vis (τ 1 ) and m c = m vis (τ 2 ), where m vis denotes the invariant mass of the measured leptons or τ h .Under this notation, m T2 is defined as: where the constraint in the minimization is over the measured lepton momenta and the missing transverse momentum, i.e. p T Σ = p T vis (τ 1 ) + p T vis (τ 2 ) + p miss T .In Eq. ( 3), the transverse mass m T is defined as and the "transverse energy" e of a particle of transverse momentum p T and mass m is defined as e = m 2 + p 2 T . ( We use the implementation in Ref. [65] to perform the minimization of Eq. ( 3).
The m T2 variable has a large discriminating power between the HH signal and the tt background, as it is bounded above by the top quark mass m t for the irreducible background process tt → bb WW → bb τν τ τν τ , while it can assume larger values for the HH signal where the tau and the b jet do not originate from the same parent particle.Detector resolution effects and other decay modes of the tt system (e.g.jets from the W boson misidentified as τ h ) result in an extension of the tail of the m T2 distribution in tt events beyond the m t value.

Background estimation
The main background sources that contaminate the signal region are tt production, Z/γ * → production and QCD multijet events.
The backgrounds from tt, single top, single Higgs boson, W boson in association with jets, and diboson processes are estimated from simulation, as described in Section 3.
The Z/γ * → background contribution is estimated using the simulation, where the LO modelling of jet emission in the Z/γ * process is known to be imperfect [66].Therefore, correction factors are calculated using events containing two isolated, opposite-sign muons compatible with the Z → µµ decay in association with two jets that satisfy similar invariant mass criteria as in the signal region.This Z+2 jets sample is divided into three control regions according to the number of b-tagged jets (0, 1, and 2) and three correction factors are derived for the Z/γ * production in association with 0, 1, or ≥2 generator level jets initiated by b quarks, and applied in the signal regions.
The multijet background is determined from data in a jet-enriched region defined by requiring that the two selected τ lepton candidates have the same electric charge.The yield is obtained from this same-sign (SS) region, where all the other selections are applied as in the signal region.The events in this region are scaled by the ratio of opposite-sign (OS) to SS event yields obtained in a multijet-enriched region with inverted τ lepton isolation.The contributions of other backgrounds, based on predictions from simulated samples, are subtracted in the OS and SS regions.The shape of the multijet background is estimated using the events in an SS region with relaxed τ lepton isolation, after subtracting the other background contributions.

Systematic uncertainties
The effects of an imperfect knowledge of the detector response, discrepancies between simulation and data, and limited knowledge of the background and signal processes are accounted for in the analysis as systematic uncertainties.They are separately treated as "normalization" uncertainties or "shape" uncertainties; the first affect the number of expected events in the signal region, while the second affect their distributions.

Normalization uncertainties
The following normalization uncertainties are considered: • The integrated luminosity is known with an uncertainty of 2.5% [67].This value is obtained from dedicated Van der Meer scans and the stability of detector response during the data taking.The uncertainty is applied to the signal and to tt, W+jets, single top quark, single Higgs boson, and diboson backgrounds, but it is not applied to the multijet and Z+jets backgrounds because they are estimated or corrected from data.• Electron, muon, and τ h lepton trigger, reconstruction and identification efficiencies are measured using Z → ee, Z → µµ, and Z → ττ → τ h ν τ µν µ ν τ events collected at √ s = 13 TeV.The corresponding uncertainties are considered as uncorrelated among the final states and are about 3% for electrons, 2% for muons, and 6% for τ leptons.
• The uncertainty in the knowledge of the τ h energy scale is about 3% for each τ h candidate [53], and its impact on the overall normalization ranges from 3 to 10% depending on the process being considered.This effect is fully correlated with a corresponding shape uncertainty in the distribution of m T2 and m KinFit HH .
• Uncertainties arising from the imperfect knowledge of the jet and b jet measured energy [50] have an impact of about 2% for the signal processes and 4% for the backgrounds.• Uncertainties in the b tagging efficiency in the simulation are evaluated as functions of jet p T and η [57] and result in an average value of 2 to 6% for the samples with genuine b jets in the final state.• For the tt process, the uncertainty in the normalization of the cross section is +4.8%/−5.5%.For the W+jets, single top quark, diboson, and single Higgs backgrounds, uncertainties range from 1 to 10%.• The uncertainties in the three correction factors derived in the control regions with 0, 1, and 2 b-tagged jets for the Z/γ * → background are propagated from the control regions to the signal region, taking into account the correlation between them, and amount to an uncertainty in the range 0.1-2.5% • The uncertainty in the multijet background normalization is estimated by propagating the statistical uncertainties in the number of events used for its determination in the region with the sign requirement inverted, as described in Section 6, and ranges between 5 and 30% depending on the final state and category.Additional sources of systematic uncertainties were found to be negligible with respect to the statistical component given the number of events in the signal and control regions.• The uncertainties in the signal cross section arising from scale variations result in an uncertainty in its normalization of +4.3%/−6.0%while effects from other theoretical uncertainties such as uncertainties on α s , PDFs and finite top quark mass effects at NNLO amount to a further 5.9% uncertainty.
The systematic uncertainties are summarized in Table 1.

Shape uncertainties
The following shape uncertainties are considered: • The shape uncertainty affecting the kinematic distribution in the simulation of the  [68], and has an impact smaller than 1% on the sensitivity of the measurement.• Uncertainties due to the limited number of simulated events or due to the statistical fluctuations of events in the multijet control region are taken into account.These uncertainties are uncorrelated across bins in the individual template shapes and their inclusion has an impact on the sensitivity smaller than 7%.• Uncertainties due to the τ h and jet energy scales are taken into account and are fully correlated with the associated normalization uncertainties.Uncertainties in the energy scales for other objects have negligible impacts on the simulated event distributions and are not taken into account.

Results
Figures 2, 3, and 4 show the distributions of the m KinFit HH and m T2 variables in the τ µ τ h , τ e τ h , and τ h τ h final states, respectively.The expected signature of resonant HH production is a localized excess in the m KinFit HH distribution, while an enhancement in the tails of the m T2 distribution would reveal the presence of nonresonant HH production.A binned maximum likelihood fit is performed simultaneously in the signal regions defined in this search for the three final states considered.The systematic uncertainties discussed previously in Section 7 are introduced as nuisance parameters in the maximum likelihood fit.In the absence of evidence for a signal, we set 95% CL upper limits on the cross section for Higgs boson pair production using the asymptotic modified frequentist method (asymptotic CL s ) [69,70].
For the resonant production mode, limits are set as a function of the mass of the resonance m S under the hypothesis that its intrinsic width is negligible compared to the experimental resolution.The observed and expected 95% CL limits are shown in Fig. 5, upper panel.The figure also shows the expectation for radion production, a spin-0 state predicted in WED models, for the parameters Λ R = 3 TeV (mass scale) and kl = 35 (size of the extra dimension), and assuming the absence of mixing with the Higgs boson.The corresponding cross section and branching fractions are taken from [71].These model-independent limits are also interpreted in the hMSSM scenario [72,73], that is a parametrization of the MSSM that considers the observed 125 GeV Higgs boson as the lighter scalar predicted from the model (usually denoted as h in the context of the model), while the resonance of mass m S represents the heavier CP-even scalar (usually denoted as H in the context of the model).Excluded regions as a function of the         For the nonresonant production mode, including the theoretical uncertainties, the observed 95% CL upper limit on the HH production cross section times branching fraction amounts to 75.4 fb while the expected 95% CL upper limit amounts to 61.0 fb.These values correspond to about 30 and 25 times the SM prediction, respectively.Limits are set for different hypotheses of anomalous self-coupling and top quark coupling of the Higgs boson.The signal kinematics depend on the ratio of the two couplings and 95% CL upper limits are set as a function of k λ /k t , assuming the other BSM couplings to be zero.The result is shown in Fig. 6, upper panel, and the exclusion is compared with the theoretical prediction for the cross section for k t = 1 and k t = 2.The sensitivity varies as a function of k λ and k t because of the corresponding changes in the signal m T2 distribution.These upper limits are used to set constraints on anomalous k λ and k t couplings as shown in Fig. 6, lower panel, where the c 2 , c 2g , and c g couplings are assumed to be equal to zero.The branching fractions for the decays of the Higgs boson into a bb and ττ pair are assumed to be those predicted by the SM for all the values of k λ and k t tested.

Summary
A search for resonant and nonresonant Higgs boson pair (HH) production in the bbττ final state is presented.This search uses a data sample collected in proton-proton collisions at √ s = 13 TeV that corresponds to an integrated luminosity of 35.9 fb −1 .The three most sensitive decay channels of the τ lepton pair, requiring the decay of one or both τ leptons into final-state hadrons and a neutrino, are used.The results are found to be statistically compatible with the expected standard model (SM) background contribution, and upper limits at the 95% confidence level are set on the HH production cross sections.
For the resonant production mechanism, upper exclusion limits at 95% confidence level (CL) are obtained for the production of a narrow resonance of mass m S ranging from 250 to 900 GeV.These model-independent results are interpreted in the context of the hMSSM scenario, where a region in the parameter space corresponding to values of m A between 230 and 360 GeV and tan β 2 is excluded at 95% CL.
For the nonresonant production mechanism, the theoretical framework of an effective Lagrangian is used to parametrize the cross section as a function of anomalous couplings of the Higgs boson.Upper limits at 95% CL on the HH cross section are obtained as a function of k λ = λ HHH /λ SM HHH and k t = y t /y SM t .The observed 95% CL upper limit corresponds to approximately 30 times the theoretical prediction for the SM cross section, and the expected limit is about 25 times the SM prediction.This is the highest sensitivity achieved so far for SM HH production at the LHC.

IG. 2 :
Feyman diagrams contributing to double Higgs production via gluon fusion (an additional ontribution comes from the crossing of the box diagram).The last diagram on the first line ontains the tthh coupling, while those in the second line involve contact interactions between the iggs and the gluons denoted with a cross.
FIG.1: Cartoon of the region in the plane (g ⇤ , /g ⇤ ), defined by by an analysis including only dimension-6 operators (in white) description is possible in the gray area ( < g min ), while ex (g min < < p g ⇤ g min ) requires including the dimension-8 opera

FIG. 2 :
FIG.2: Feyman diagrams contributing to double Higgs product contribution comes from the crossing of the box diagram).T contains the tthh coupling, while those in the second line involv Higgs and the gluons denoted with a cross.

IG. 2 :
Feyman diagrams contributing to double Higgs production via gluon fusion (an additional ontribution comes from the crossing of the box diagram).The last diagram on the first line ontains the tthh coupling, while those in the second line involve contact interactions between the iggs and the gluons denoted with a cross.

Figure 2 :
Figure 2: Distributions of the events observed in the signal regions of the τ µ τ h final state.The first, second, and third rows show the resolved 1b1j, 2b, and boosted regions, respectively.Panels in the right column show the distribution of the m T2 variable, while the other panels show the distribution of the m KinFit HH variable, separated in the low-mass (LM, left panels) and high-mass (HM, central panels) regions for the resolved event categories.Data are represented by points with error bars and expected signal contributions are represented by the solid (BSM HH signals) and dashed (SM nonresonant HH signal) lines.Expected background contributions (shaded histograms) and associated systematic uncertainties (dashed areas) are shown as obtained after the maximum likelihood fit to the data under the background-only hypothesis.The background histograms are stacked while the signal histograms are not stacked.

Figure 3 :
Figure 3: Distributions of the events observed in the signal regions of the τ e τ h final state.The first, second, and third rows show the resolved 1b1j, 2b, and boosted regions, respectively.Panels in the right column show the distribution of the m T2 variable, while the other panels show the distribution of the m KinFit HH variable, separated in the low-mass (LM, left panels) and high-mass (HM, central panels) regions for the resolved event categories.Data are represented by points with error bars and expected signal contributions are represented by the solid (BSM HH signals) and dashed (SM nonresonant HH signal) lines.Expected background contributions (shaded histograms) and associated systematic uncertainties (dashed areas) are shown as obtained after the maximum likelihood fit to the data under the background-only hypothesis.The background histograms are stacked while the signal histograms are not stacked.

Figure 4 :
Figure 4: Distributions of the events observed in the signal regions of the τ h τ h final state.The first, second, and third rows show the resolved 1b1j, 2b, and boosted regions, respectively.Panels in the left column show the distribution of the m KinFit HH variable and panels in the right column show the distribution of the m T2 variable.Data are represented by points with error bars and expected signal contributions are represented by the solid (BSM HH signals) and dashed (SM nonresonant HH signal) lines.Expected background contributions (shaded histograms) and associated systematic uncertainties (dashed areas) are shown as obtained after the maximum likelihood fit to the data under the background-only hypothesis.The background histograms are stacked while the signal histograms are not stacked.

m
A and tan β parameters, representing respectively the mass of the CP-odd scalar and the ratio of the vacuum expectation values of the two Higgs doublets of the model, are shown in Fig.5, lower panel.The minimum of the sensitivity around m S = 270 GeV results in the presence of two separate expected excluded regions in this interpretation.

Figure 5 :Figure 6 :
Figure5: (upper) Observed and expected 95% CL upper limits on cross section times branching fraction as a function of the mass of the resonance m S under the hypothesis that its intrinsic width is negligible with respect to the experimental resolution.The inner (green) band and the outer (yellow) band indicate the regions containing 68 and 95%, respectively, of the distribution of limits expected under the background-only hypothesis.The red line denotes the expectation for the production of a radion, a spin-0 state predicted in WED models, for the parameters Λ R = 3 TeV (mass scale) and kl = 35 (size of the extra dimension), assuming the absence of mixing with the Higgs boson.(lower) Interpretation of the exclusion limit in the context of the hMSSM model, parametrized as a function of the tan β and m A parameters.In this model, the CP-even lighter scalar is assumed to be the observed 125 GeV Higgs boson and is denoted as h, while the CP-even heavier scalar is denoted as H and the CP-odd scalar is denoted as A. The dotted lines indicate trajectories in the plane corresponding to equal values of the mass of the CP-even heavier scalar of the model, m H .

5 Signal regions and discriminating observables signal
and background: ∆ϕ(H bb , H ττ ), ∆ϕ(H ττ , p miss ∆R(b, b) p T (H bb ), ∆R( , τ h ) p T (H ττ ), m T ( ), m T (τ h ), and ∆ϕ( , p miss Here refers to the selected muon or electron, H bb and H ττ denote the H boson candidates reconstructed from the two jets and the two τ leptons, respectively, and m T ( ) = √ (p T + p miss T ) 2 − ( p T + p miss T T), ∆ϕ(H bb , p miss T ),T).

Table 1 :
Systematic uncertainties affecting the normalization of the different processes.estimated by varying the top quark p T distribution according to the uncertainties in differential p T measurements described in Ref.
Individuals have received support from the Marie-Curie programme and the European Research Council and Horizon 2020 Grant, contract No. 675440 (European Union); the Leventis Foundation; the A. P. Sloan Foundation; the Alexander von Humboldt Foundation; the Belgian Federal Science Policy Office; the Fonds pour la Formation à la Recherche dans l'Industrie et dans l'Agriculture (FRIA-Belgium); the Agentschap voor Innovatie door Wetenschap en Technologie (IWT-Belgium); the Ministry of Education, Youth and Sports (MEYS) of the Czech Republic; the Council of Science and Industrial Research, India; the HOMING PLUS programme of the Foundation for Polish Science, cofinanced from European Union, Regional Development Fund, the Mobility Plus programme of the Ministry of Science and Higher Education, the National Science Center (Poland), contracts Harmonia 2014/14/M/ST2/00428, Opus 2014/13/B/ST2/02543, 2014/15/B/ST2/03998, and 2015/19/B/ST2/02861, Sonata-bis 2012/07/E/ST2/01406; the National Priorities Research Program by Qatar National Research Fund; the Programa Clarín-COFUND del Principado de Asturias; the Thalis and Aristeia programmes cofinanced by EU-ESF and the Greek NSRF; the Rachadapisek Sompot Fund for Postdoctoral Fellowship, Chulalongkorn University and the Chulalongkorn Academic into Its 2nd Century Project Advancement Project (Thailand); and the Welch Foundation, contract C-1845.