Searches for Resonant Scalar Boson Pair Production Using Run 2 LHC Proton-Proton Collision Data

The discovery of the Higgs boson in 2012 provided confirmation of spontaneous electroweak symmetry breaking as the mechanism by which fundamental particles gain mass and thus completed the Standard Model of particle physics. Additionally, it opened a new approach to searching for potential new particles. Many beyond the Standard Model theories predict new heavy particles that couple to the Higgs boson, leading to a resonant production mode of Higgs boson pairs. Other theories extend the Higgs sector by introducing additional scalar bosons that differ from the observed Higgs boson only by mass. The ATLAS and CMS Collaborations have searched for evidence of such processes using s=13 TeV Run 2 proton-proton collision data at the Large Hadron Collider. This review article summarizes the latest experimental results from searches for resonant production of pairs of Higgs bosons or additional Higgs-like scalar bosons at ATLAS and CMS.


Introduction
The 2012 discovery of the Higgs boson (H) [1,2] provided the final missing piece of the Standard Model (SM) of particle physics by confirming the mechanism of spontaneous electroweak symmetry breaking to generate masses for fundamental particles. As with the discovery of all new particles, the discovery of the Higgs boson has provided an important new means of searching for evidence of new Beyond the Standard Model (BSM) physics effects or additional heavy particles. Many BSM theories predict the existence of new heavy particles that couple strongly to H, leading to the potential resonant production of Higgs boson pairs (HH). The observation of such a resonant production would appear as an excess in HH production over the SM prediction localized in the HH invariant mass (m HH ) spectrum.
An overview of searches for resonant production of scalar bosons by the ATLAS and CMS Collaborations performed using √ s = 13 TeV Run 2 proton-proton (pp) collision data at the Large Hadron Collider (LHC) is presented here. The searches use either a partial Run 2 dataset or the full Run 2 dataset. A comprehensive review of the status of theoretical and experimental efforts related to searches for HH production was previously presented in Reference [3]. The searches presented in this review include several updated results since the comprehensive review.
In the descriptions of all searches provided in this review, the term "lepton" refers to an electron (e) or a muon (µ), unless otherwise specified. Furthermore, notation varies between analyses. Unless specified otherwise in an analysis-specific context, H refers to a SM Higgs boson or the observed Higgs boson with a mass of approximately 125 GeV, S refers to an additional scalar that has mass-dependent couplings similar to those of the Higgs boson, and X refers to a heavy scalar.
The searches in this review all use simulated signal samples. Except where specified, these simulated samples are generated assuming no interference with SM processes and assuming 100% branching fractions into SM decay products. Scalar resonances simulated with these assumptions are referred to as "generic scalars" in this text. Resonances are typically generated with a width that is smaller than detector resolution, and interference effects from SM HH production are neglected.

Theoretical Models
A brief overview of several theoretical models that predict the resonant production of scalar bosons is presented here. It is not meant to be a comprehensive list of models, but covers the most relevant models for the experimental searches that are discussed. More details on each model presented here can be found in the respective references.
In all UV complete renormalizable models, couplings of charge-parity-even (CP-even) scalars (H i ) to vector bosons (V) are modified from those of the SM Higgs boson with the sum rule imposed by unitarity constraints [4]. The coupling of the observed Higgs boson to vector bosons is consistent with the SM predictions, resulting in strongly suppressed couplings for the other scalars predicted by such models.

Resonant HH Production
There are many BSM theories that predict new heavy particles that can potentially be created at the LHC and are expected to decay to HH. These new particles can either have a spin of 0 or 2.

Singlet Models
The simplest extension of the SM Higgs sector is the addition of a real gauge singlet scalar field [5][6][7][8][9]. In general, the singlet field mixes with the SM Higgs boson, resulting in couplings to SM particles. After electroweak symmetry breaking, such extensions result in two scalar bosons (H and X). Assuming m X > 2m H where m H ≈ 125 GeV, the mixing can allow for X → HH decays, resulting in the resonant production of a pair of on-shell H. Mixing between singlets and the SM doublet can result in significant interference effects in the resonant and SM HH production modes [10]. Singlet benchmark models are discussed in detail in References [3,[11][12][13]].

2HDM
One of the most common BSM extensions to the Higgs sector is the addition of second Higgs doublet, leading to a class of models referred to as Two-Higgs-Doublet Models (2HDMs) [14,15]. There are several types of 2HDMs, classified by the couplings of each Higgs doublet to fermions. In Type-I 2HDMs, one doublet couples to up-and down-type quarks as well as charged leptons while the other doublet does not couple to fermions. In Type-II 2HDMs, one doublet couples to up-type quarks, and the other couples to downtype quarks and charged leptons. In lepton-specific 2HDMs, one doublet couples to all quarks and the other couples to charged leptons. In flipped-type 2HDMs, one doublet couples to up-type quarks and charged leptons, and the other couples to down-type quarks. There are additionally Type-III 2HDMs in which both doublets couple to all fermions and may or may not allow tree-level flavor changing neutral currents.
All 2HDM extensions of the SM result in a total of five Higgs bosons. If there is no CP violation, these consist of two neutral CP-even scalars (H and X), a neutral CP-odd pseudoscalar (A), and two charged bosons (H ± ). Typically, H is treated as the observed Higgs boson with m H ≈ 125 GeV, and X is more massive and is often used as a benchmark spin-0 resonance that decays to HH. If there is CP violation, the three neutral states mix to form mass eigenstates with indefinite CP quantum numbers, one of which is identified as the observed Higgs boson [16]. Multiple finals states are used in the searches presented here to make use of the 159 wide variety of Higgs boson decay modes with large branching ratios or clean signatures. 160 The HH branching ratios are determined from the known branching ratios of the SM 161 Higgs boson [13]. A summary of the branching ratios of the most important HH decay 162 Figure 1. Feynman diagrams depicting the ggF production mode for (a) a scalar resonance and (b) a Kaluza-Klein graviton that subsequently decay into a pair of Higgs bosons [42]. Multiple finals states are used in the searches presented here to make use of the wide variety of Higgs boson decay modes with large branching ratios or clean signatures. The HH branching ratios are determined from the known branching ratios of the SM Higgs boson [13]. A summary of the branching ratios of the most important HH decay modes is shown in Figure 3. Each final state has advantages and disadvantages that determine its sensitivity across kinematic regimes. For a universally optimal search, a combination of many final states is required.  Branching ratios for the most important HH decay modes assuming SM couplings, calculated at NLO, from Reference [13]. [Image credit: K. Leney].

Searches for HH Production in the bbbb Channel
The Standard Model Higgs boson decays predominantly to bb with a branching ratio of BR(H → bb) ≈ 58%. As a result, HH → bbbb is the dominant HH decay mode with a branching ratio of BR(HH → bbbb) ≈ 33%. This large branching ratio makes bbbb an important search channel, especially for resonance masses above 1 TeV. The purely hadronic signature results in a large multijet background that is mitigated through the use of b-tagging methods to identify hadronic jets consistent with the decay of B-hadrons.

ATLAS bbbb Searches
The ATLAS Collaboration has conducted two searches for resonant HH production in the bbbb final state using the full Run 2 dataset. This corresponds to 126-139 fb −1 of data collected between 2015 and 2018. The analyses are designed to search for resonances in the ggF and VBF production modes.

ATLAS ggF bbbb Search
The first analysis [42] is a search for ggF resonant HH production, and the results are interpreted in the context of spin-0 and spin-2 benchmark models. The heavy CP-even neutral scalar in a 2HDM is used as the spin-0 benchmark model and is generated with a narrow width that is much smaller than the detector resolution. The spin-2 benchmark model is a Kaluza-Klein graviton with κ/M Pl = 1 with a generated width ranging from 3% to 13% of the resonance mass.
Both resolved and boosted topologies are exploited to set limits on the two models in the mass range from 251 GeV to 3 TeV. The resolved channel is used to search for resonances in the mass range from 251 GeV to 1.5 TeV, and the boosted channel is used to search for resonances in the mass range from 900 GeV to 3 TeV. Both channels are statistically combined in the overlapping mass range.
In the resolved channel, pairs of H candidates are reconstructed from four anti-k t R = 0.4 jets using a boosted decision tree (BDT). Events are split into 4b and 2b categories if they contain at least four b-tagged or exactly two b-tagged jets. In the 4b category, the four p T -leading b-tagged jets are used, and in the 2b category, the two b-tagged jets and the two p T -leading non-b-tagged jets are used. Events are selected using a combination of triggers requiring one or two high-E T b-tagged jets plus additional jets from 126 fb −1 of data collected between 2016 and 2018. An upper bound cut is applied on the pseudorapidity separation between the two H candidates to reduce the quantum chromodynamics (QCD) multijet background. Additionally, a cut is applied to reject any events that are kinematically consistent with the hadronic decay of a top quark.
Three regions are defined in the m(H 1 ) − m(H 2 ) plane, where m(H 1 ) and m(H 2 ) denote the masses of the reconstructed Higgs boson candidates ordered in p T . The signal region (SR), control region (CR), and validation region (VR) are shown in Figure 4, superimposed over data in the 2b category.  Following the event selection, the "corrected m(HH)" is constructed as the final discriminating variable. This is performed by rescaling the four momenta of the two H candidates to obtain m(H 1 ) = m(H 2 ) = 125 GeV.
The background after the selection cuts are applied is approximately 95% composed of QCD multijet events and is estimated using a data-driven method. Other background processes are ignored. A neutral network (NN) is used to estimate the background in the 4b SR from data in the 2b SR. The corrected m(HH) distribution in the 4b SR is shown in Figure 5 for the background estimate, data, and several spin-0 signal hypotheses. Data Background Background Figure 5. Distribution of the corrected m(HH) for data, the background estimate, and select spin-0 signal hypotheses in the resolved 4b SR after the background-only fit. The background distribution is shown with its total post-fit uncertainty. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The signal hypotheses are normalized to their expected cross-section exclusion limits. The bottom panel shows the relative difference between the data and background estimate [42].
In the boosted channel, two anti-k t R = 1.0 (large-R) jets are used to reconstruct the two H candidates. The channel uses 139 fb −1 of data collected between 2015 and 2018. Events are selected using a combination of triggers that require at least one high-E T large-R jet.
Events are required to contain at least two large-R jets with a mass of m J > 50 GeV and have at least one ghost-associated [43] variable radius track jet. The two p T -leading large-R jets that satisfy the mass and track jet requirements are selected as the two H candidates. An upper bound cut is applied on the pseudorapidity separation between the two H candidates to reduce the QCD multijet background.
The events are split into three signal-enriched (high-tag) categories based on the number of b-tagged track jets associated with each of the H candidates. The 4b category requires each H candidate to have at least two associated b-tagged track jets. The 3b category requires one H candidate to have at least two associated b-tagged track jets and the other H candidate to have exactly one associated track jet. The 2b category requires each H candidate to have exactly one associated b-tagged track jet. Each high-tag category has an associated "low-tag" category that is used to estimate the background. Sketches of H candidates in the high-tag and low-tag categories are shown in Figure 6. 2b-2f 2b-1f 1b-1f Figure 6. Sketches of the three signal-enriched and three low-tag categories. The teal cones represents the H candidate large-R jets, the yellow cones represent associated b-tagged track jets, and white cones represent associated untagged track jets [42].
Similar to the resolved channel, the SR, CR, and VR are defined in the m(H 1 ) − m(H 2 ) plane. The three regions are shown in Figure 7 superimposed over data from events from the 3b "low-tag" category, denoted as 2b-1 f . The background after the selection cuts are applied is dominated by QCD multijet events with a tt contribution ranging from 10% to 30% of the background. Other processes contribute ≤ 1% of the background and therefore are neglected. The tt background is estimated using Monte Carlo (MC) simulation with data-driven corrections. The QCD multijet background is estimated using a data-driven technique. A normalization fit in the CR data is used to simultaneously derive normalization factors for the QCD multijet and tt backgrounds in each b-tagging category. A kinematic reweighting procedure based on iterative splines is used to estimate the background in the SRs from the "low-tag" categories. The m(HH) distributions in the 2b and 4b SRs are shown in Figure 8 for the background estimate, data, and several spin-0 signal hypotheses. The background distribution is shown with its total post-fit uncertainty. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. Select spin-0 signal hypotheses are overlaid, normalized to a cross-section equal to the expected cross-section exclusion. The bottom panel shows the relative difference between the data and background model. [42] prediction for the bulk RS model with κ/M Pl = 1 is also shown. The model is excluded 263 for masses between 298 GeV and 1440 GeV. The second analysis [45] searches for resonant HH production in the VBF produc-266 tion mode and the results are interpreted in the context of spin-0 resonance benchmark 267 models. A Type II 2HDM [15] is used as the broad resonance with a width of 10-20% of 268 the resonance mass and is obtained by setting the ratio of the vacuum expectation values 269 for the two Higgs doublet to tan(β) = 2.0 and sin(β − α) = 0.6 where α is the mixing 270 angle between the two CP-even Higgs bosons. It should be noted that this parameter 271 choice has been excluded by a combination of Higgs boson measurements [46]. A generic 272 scalar with a fixed width of 4 MeV is used as the narrow resonance.   Figure 8. Distributions of m(HH) for data, the background estimates, and example signal hypotheses in the boosted (a) 2b and (b) 4b SRs after the background-only fit. The backgrounds distribution are shown with their total post-fit uncertainty. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The signal hypotheses are normalized to their expected cross-section exclusion limits. The bottom panels show the relative differences between the data and background estimates [42].
Systematic uncertainties are evaluated for both channels. The dominant source of uncertainty for both is the data-driven background estimate. Several components are evaluated to cover various aspects of the methods. Additionally, experimental uncertainties are evaluated for the detector modeling, especially of hadronic jets and the b-tagging methods. Finally, theoretical systematic uncertainties are evaluated on the renormalization and factorization scales and parton showering for the signal. These are also varied for the tt background simulation along with the POWHEG damping parameter, which regulates radiation at high transverse momentum.
A profile likelihood ratio test statistic is used to test the production cross-section times branching ratio for the signal models. The likelihood fit is conducted in bins of corrected m(HH) for the resolved channel and m(HH) for the boosted channel. Data are fit separately in the resolved and boosted channels. All categories are fit simultaneously in the boosted channel. In the boosted channel, the 2b category is only used for resonance masses of at least 2 TeV. No significant excess over the background prediction is observed.
Upper limits on the production cross-section times branching ratio are set at the 95% confidence level using the asymptotic CL s method [44]. The limits from both channels are combined statistically in the resonance mass range of 900 GeV to 1.5 TeV. The combined expected and observed limits are shown in Figure 9. The theoretical prediction for the bulk RS model with κ/M Pl = 1 is also shown. The model is excluded for masses between 298 and 1440 GeV. ATLAS Preliminary s = 13 TeV, 126 − 139 fb −1 Spin-0 Observed limit (95% CL) Expected limit (95% CL) Expected limit ± 1σ Expected limit ± 2σ Resolved expected limit Boosted expected limit ATLAS Preliminary s = 13 TeV, 126 − 139 fb −1 Spin-2 Observed limit (95% CL) Expected limit (95% CL) Expected limit ± 1σ Expected limit ± 2σ Resolved expected limit Boosted expected limit RS Graviton, k/M Pl = 1 (b) Spin-2 model Figure 9. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant ggF HH production in the bb¯bb¯ final state using the full Run-2 ATLAS dataset. Limits are set on benchmark (a) spin-0 and (b) spin-2 models. Expected limits are shown individually for the resolved and boosted channels, as well as combined expected and observed limits. The ±1σ and ±2σ uncertainties on the expected limits are shown with the colored bands. The theoretical prediction for the RS model with κ/M Pl = 1 is shown in (b). [42] Figure 9. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant ggF HH production in the bbbb final state using the full Run 2 ATLAS dataset. Limits are set on benchmark (a,(b). Expected limits are shown individually for the resolved and boosted channels, as well as combined expected and observed limits. The ±1σ and ±2σ uncertainties on the expected limits are shown with the colored bands. The theoretical prediction for th e RS model with κ/M Pl = 1 is shown in (b) [42].

ATLAS VBF bbbb Search
The second analysis [45] is a search for VBF resonant HH production, and the results are interpreted in the context of spin-0 resonance benchmark models. A Type II 2HDM [15] with a width of 10-20% of the resonance mass is used as the broad resonance benchmark mode. It is obtained by setting the ratio of the vacuum expectation values for the two Higgs doublet to tan(β) = 2.0 and sin(β − α) = 0.6, where α is the mixing angle between the two CP-even Higgs bosons. It should be noted that this parameter choice has been excluded by a combination of Higgs boson measurements [46]. A generic scalar with a fixed width of 4 MeV is used as the narrow resonance benchmark model. The analysis uses resolved reconstruction techniques and sets limits on both models in the mass range of 260 GeV to 1 TeV. Pairs of H candidates are reconstructed from b-tagged central anti-k t R = 0.4 jets, and forward jets are used to select events with the VBF topology.
Events are selected using a combination of b-jet triggers requiring at least one or two b-tagged jets from 126 fb −1 of √ s = 13 TeV data collected between 2016 and 2018. Events are required to have exactly four b-tagged central jets and at least two forward jets. The two p T -leading forward jets with an opposite sign of η, referred toas VBF jets, are required to have a large pseudorapidity separation and have an invariant mass greater than 1 TeV. The four central b-tagged jets are paired into the three possible combinations to construct the two H candidates. To be considered, each pairing is required to satisfy requirements on the ∆R between each pair of jets. If there is more than one pairing that satisfies the ∆R requirements, the pairing that results in invariant masses that are the most consistent with two H is used.
Background events are suppressed with requirements designed to select signal-like events. Requirements are placed on the p T of the vector sum of the jets, the pseudorapidity between the H candidates, and a variable (X Wt ) designed to reject events that are kinematically consistent with the hadronic decay of a top quark.
Similar to the ggF analysis, three regions are defined in the m lead 2b − m subl 2b plane, where m lead 2b and m subl 2b are the masses p T -leading and p T -subleading H candidates, respectively. The SR, VR, and sideband region (SB) are shown in Figure 10 superimposed over the multijet background estimate.  Following the event selection, the background consists of 95% QCD multijet events and 5% tt events. The SM HH produced via ggF is also considered as a small source of background events. Other background sources are neglected. The tt and SM ggF HH backgrounds are estimated using MC. The QCD multijet background is estimated using a data-driven technique.
The QCD multijet background is estimated using data in a region identical to the SR but with exactly two b-tagged central jets and at least two untagged central jets. Corrections are Symmetry 2022, 14, 260 12 of 66 applied to account for kinematic differences due to the different b-tagged jet multiplicities, following the procedure described in Reference [47].
The QCD multijet and all-hadronic tt backgrounds are normalized with a simultaneous fit of the X Wt distributions to data in the SB. The non-all-hadronictt and SM ggF HH backgrounds are normalized to their SM predictions. The post-fit m 4b distribution in the SR is shown in Figure 11 for the background estimate, data, and signal hypotheses for a narrow spin-0 resonance and SM VBF HH.  . Distribution of the m 4b variable for data, the background estimate, a spin-0 signal hypothesis, and a SM HH signal hypothesis in the SR after the background-only fit. The background distribution is shown with its total post-fit uncertainty. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The signal hypotheses are normalized to their expected cross-section exclusion limits. The bottom panel shows the ratio of the data to the background estimate [45].
Systematic uncertainties are evaluated on the background estimate, theoretical predictions, and detector modeling. The QCD multijet background estimate is the source of the dominant systematic uncertainties. The experimental modeling of the jet energy scale and resolution also contribute significantly. The statistical uncertainty is larger than the combined systematic uncertainty.
A profile-likelihood ratio test statistic is used to test the production cross-section times branching ratio for the signal models. No significant excess over the background prediction is observed.
Upper limits are set at 95% confidence level using the asymptotic CL s method with m 4b as the final discriminant. Figure 12 shows the limits for the narrow and broad spin-0 models in the resonance mass range from 260 GeV to 1 TeV. Observed limit (95% CL) Expected limit (95% CL) Broad spin-0 model Figure 12. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant VBF HH production in the bb¯bb¯ final state using the full Run-2 ATLAS dataset. Limits are set on benchmark (a) narrow and (b) broad spin-0 models. The ±1σ and ±2σ uncertainties on the expected limits are shown with the colored bands. [45] Figure 12. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant VBF HH production in the bbbb final state using the full Run 2 ATLAS dataset. Limits are set on benchmark (a,b). The ±1σ and ±2σ uncertainties on the expected limits are shown with the colored bands [45].

CMS bbbb Searches
The CMS Collaboration has conducted two searches for ggF resonant HH production in the bbbb final state using Run 2 LHC data. The analyses are designed to search for signal hypotheses with low masses and resolved topologies and with high masses and boosted topologies.

CMS Resolved bbbb Search
The first analysis [48] is a search for ggF resonant HH production with a resolved topology using a partial Run 2 dataset corresponding to 35.9 fb −1 of data collected in 2016.
The search is interpreted in the context of narrow spin-0 and spin-2 benchmark models. A bulk radion with Λ R = 3 TeV and κl = 35 is used as the spin-0 benchmark model and a Kaluza-Klein graviton with κ/M Pl = 0.5 and κl = 35 [49] is used as the spin-2 benchmark model. Both benchmark models are generated with a width of 1 MeV. Limits are set on both benchmark models in the mass range from 260 GeV to 1.2 TeV. The analysis is split into the low-mass region (LMR) for resonance masses less than 580 GeV and the medium-mass region (MMR) for resonance masses greater than 580 GeV.
Events are selected with a combination of two triggers that require at least four anti-k t R = 0.4 jets, at least three of which are required to be b-tagged. After the trigger selection, events are required to contain at least four b-tagged jets. The selected jets are paired to form A multivariate regression technique is used to correct the b-tagged jets' p T , following the procedures in References [50,51]. This improves the mass resolution of the reconstructed H candidates from 10-13% to 6-12% as well as the mass scale. An additional correction is applied to the b-tagged jet momenta to constrain the H candidate masses to be 125 GeV. The impact of the regression and constraint on the resonance mass, reconstructed from the two H candidates, for spin-2 resonances is shown in Figure 13.  Figure 13. The reconstructed heavy resonance mass distribution for 450, 750, and 1000 GeV spin-2 mass hypotheses. The distributions are shown without any corrections, with the kinematic constraint on m H applied and with the regression correcting b-tagged jet energy additionally applied. The distributions are all normalized to have the same integral [48].
Two regions are defined in the m H 1 − m H 2 plane based on R. The SR and SB are shown in Figure 14 superimposed over data from events in the MMR. The signal and background are both estimated in the SR using parametric models, fit separately in the LMR and MMR. The signal parametric models are obtained using a Gaussian function in the MMR and the sum of two Gaussian functions in the LMR fit to the simulated distribution of the reconstructed resonance mass m X . The background consists primarily of QCD multijet events with a 10-15% contribution from tt, which has an m X shape similar to that of the QCD multijet background. The parametric model used to describe the background is fit to the m X distribution of data in control regions consisting of the SB as well as events in both the SR and SB in which one of the four jets used to reconstruct m X fails the b-tagging criteria. The LMR is split into LMR1 with m X ∈ [260, 310] GeV and LMR2 with m X ∈ [310, 580] GeV. The m X distributions for data in the SR and the background model for LMR1 and MMR are shown in Figure 15.  Figure 15. Distribution of m X for data in the SR and the background estimate for the (a) LMR1 and (b) MMR. The number of degrees of freedom in each fit to the background-only hypothesis is given as n. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The lower panels show the difference between the data and fits, normalized by the uncertainty in the number of data events. [48] Systematic uncertainties are evaluated on the background estimate, the choice of 366 PDF, and modeling of the detector response.

367
The production cross-section times branching ratio for the signal models is tested 368 in bins of m X . No significant excess over the background prediction is observed. Upper 369 limits on the production cross-section times branching ratio are set on both models 370 at 95% CL using the CL s method. The limits are shown in Figure 16 Figure 15. Distribution of m X for data in the SR and the background estimate for the (a,b). The number of degrees of freedom in each fit to the background-only hypothesis is given as n. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The lower panels show the difference between data and the fits, normalized by the uncertainty in the number of data events [48].
Systematic uncertainties are evaluated on the background estimate, the choice of PDF, and modeling of the detector response.
The production cross-section times branching ratio for the signal models is tested in bins of m X . No significant excess over the background prediction is observed. The upper limits on the production cross-section times branching ratio are set on both models at 95% CL using the CL s method. The limits are shown in Figure 16

CMS Boosted bbbb Search
The second analysis [52] is a search for ggF resonant HH production with a boosted topology using the full Run 2 dataset. This corresponds to 138 fb −1 of data collected between 2015 and 2018.
The search is interpreted in the context of narrow spin-0 and spin-2 benchmark models. A bulk radion with Λ R = 3 TeV and κl = 35 is used as the spin-0 benchmark model, and a Kaluza-Klein graviton with κ/M Pl = 0.5 and κl = 35 is used as the spin-2 benchmark model. Both benchmark models are generated with a width of 1 MeV. Limits are set on both benchmark models in the mass range from 1 to 3 TeV. The analysis is split into the fully-boosted category where both H → bb are reconstructed with large-R jets and the semiresolved category where one H → bb is reconstructed with a large-R jet and the other is reconstructed with two anti-k t R = 0.4 jets.
Events are selected with a combination of triggers that require hadronic jets. The first trigger requires a high scalar sum of the p T of all anti-k t R = 0.4 jets in the event (H T ). The second trigger requires a large H T and a pair of jets with an invariant mass greater than 900 GeV and a small psuedorapidity separation. The third trigger requires an anti-k t R = 0.8 jet with a mass of at least 30 GeV.
Events in the fully-boosted category are required to have at least two anti-k t R = 0.8 jets, the p T -subleading of which is required to have a mass between 110 and 140 GeV. The DeepAK8 tagger [53] is used to identify large-R jets as coming from H → bb decays. The tagger provides either loose or tight thresholds, leading to two signal regions, one in which both large-R jets pass the tight criteria and the other in which both jets pass at least the loose criteria and no more than one jet passes the tight criteria.
Events in the semiresolved category are required to have at least one anti-k t R = 0.8 jet that passes the DeepAK8 tight threshold and at least two b-tagged anti-k t R = 0.4 jets. Events are required to have at least one pair of b-tagged R = 0.4 jets where both jets are ∆R > 0.8 from the H → bb large-R jet and are within ∆R < 1.5 from each other. If multiple pairs of b-tagged R = 0.4 jets satisfy this criteria, the pair with the largest sum of b-tagging scores is selected as the resolved H → bb candidate. This is first tried with the p T -leading large-R jet and if no pairs are found, it is tried with the p T -subleading large-R jet. The large-R jet that leads to a successful pairing is identified as the boosted H → bb candidate. If no pair is found, the event is rejected. The invariant mass of the selected pair of R = 0.4 jets is required to be 90 GeV < m jj < 140 GeV.
A modified "reduced" invariant mass of the HH system, m red , is used in place of the direct invariant mass of the selected jets. The reduced mass is defined as for the boosted category where m J J is the invariant mass of the two H → bb large-R jets and m J 1 and m J 2 are the masses of the individual large-R jets. For the semiresolved category, the reduced mass is defined as for the boosted category where m Jjj is the invariant mass of the two H → bb candidates, m J is the mass of the H → bb large-R jet, and m jj is the invariant mass of the H → bb R = 0.4 jet pair. Events in both categories are required to have m red > 750 GeV. Additionally, a maximum pseudorapidity separation requirement is imposed on the two H → bb candidates.
The background consists almost entirely of QCD multijet and tt events, both of which are estimated using data-driven methods that also make use of simulation. The QCD multijet background estimate is based on the ratio of events whose p T -leading large-R jet passes the DeepAK8 tagging criteria to events whose whose p T -leading large-R jet fails the tagging criteria. The ratio is found using simulation and corrected using data outside the mass window requirement on the p T -subleading large-R jet or R = 0.4 jet pair to predict the background as a function of m red . The tt background is estimated using simulation that is fit to data in dedicated control regions that shift the mass window requirements to contain the top quark mass peak. Slices of the m J − m red plane for data, the background estimate, and a signal hypothesis are shown in Figure 17.  Figure 17. Slices of the m J − m red distributions for data, the background estimate and the spin-0 m X = 1.5 TeV signal hypothesis projected onto the (a) m J and (b) m red axes. The hatched bands represent the total uncertainty. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The lower panels show the difference between data and the background estimate normalized by the total uncertainty. [52] are also shown. The bulk radion model is excluded in the mass range of [1, 2.6] TeV. The   Figure 17. Slices of the m J − m red distributions for data, the background estimate, and the spin-0 m X = 1.5 TeV signal hypothesis projected onto the (a) m J and (b) m red axes. The hatched bands represent the total uncertainty. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The lower panels show the difference between data and the background estimate normalized by the total uncertainty [52].
Systematic uncertainties are evaluated on the background estimate, the choice of PDF, and modeling of the detector response. The dominant uncertainties come from the tt cross-section and the QCD multijet background estimate.
A two-dimensional likelihood fit is performed in the m J − m red plane. No significant excess over the background prediction is observed. The upper limits on the production cross-section times branching ratio are set on both models at 95% CL using the CL s method. The limits are shown in Figure 18

Discussion of bbbb Searches
ATLAS and CMS both performed searches in the bbbb channel using resolved and boosted topologies. For low resonance masses, the bbbb channel presents the challenges of a large QCD multijet background and a roughly isotropic resolved signal topology that is difficult to distinguish from background. For high resonance masses, however, the large bbbb branching ratio can be efficiently exploited using boosted topologies.
In the resolved topology, the greatest difficulty comes from pairing b-tagged jets to form two H candidates. This pairing is necessary to provide kinematic features that can be used to separate signal from background. As shown in Figure 10, algorithms with strong m H requirements, such as the one used in the ATLAS VBF search, can sculpt the background distribution to peak exactly in the signal region. The use of machine learning techniques to pair the jets, such as in the ATLAS resolved ggF search, can help to reduce this sculpting while maintaining a good signal reconstruction efficiency. As shown in Figure 4, this can result in less severe sculpting of the background estimate with a broad peak partially overlapping with the signal region. The MMR in the CMS resolved search uses a ∆R requirement to determine the pairing strategy, which does not sculpt the background significantly. As shown in Figure 14, this results in a smooth data distribution, consisting primarily of QCD multijet background events.
The use of large-R jets in the boosted topology significantly simplifies the reconstruction and identification of H → bb candidates. Particle-flow large-R jets identified with the constituent-based DeepAK8 algorithm used by CMS and calorimeter-based large-R jets with b-tagged track jets used by ATLAS both provide excellent performance.Additionally, the simple reconstruction of using two large-R jets or one large-R jet plus two b-tagged small-R jets does not appreciably sculpt the background estimate as shown in Figures 7 and 17.
As shown in Figure 9, there is a significant range of resonance masses for which boosted and resolved reconstruction techniques are comparably sensitive. Ensuring orthogonality between the selection criteria and combining the two analysis techniques in the region of overlap can result in significantly stronger limits than either topology alone.

Searches for H H Production in the bbτ + τ − Channel
The bbτ + τ − channel is the most sensitive to search for resonances with masses around 500 GeV decaying to HH. The large branching ratio of H → bb combined with the relatively large branching ratio and unique signature of H → τ + τ − ensures a good signal-tobackground ratio. Both leptonically and hadronically decaying τ leptons are used to ensure nearly complete coverage of final states. The backgrounds in this channel primarily consist of events with a top quark, a Z boson, or one or more hadronic jets misidentified as an electron, muon, or hadronically decaying τ lepton.

ATLAS bbτ + τ − Searches
The ATLAS Collaboration has conducted two searches for ggF resonant HH production in the bbτ + τ − final state using the full Run 2 dataset of 139 fb −1 of data collected between 2015 and 2018. The first analysis uses resolved topologies to search for resonance masses between 251 GeV and 1.6 TeV. The second analysis uses boosted topologies to search for resonance masses between 1 and 3 TeV.
The first analysis [54] is a search for ggF resonant HH production with a resolved topology in the τ had τ had and τ lep τ had channels. The search is interpreted in the context of a generic narrow spin-0 benchmark model. The signal samples are generated with a width of 10 MeV. Limits are set on the benchmark model in the mass range from 251 GeV to 1.6 TeV.
The analysis is split into two channels based on the decay mode of the τ + τ − system. The channel in which both τ leptons decay hadronically (τ had ) is referred to as τ had τ had and the channel in which one τ lepton decays hadronically and the other decays leptonically (τ lep ) is referred to as τ lep τ had . The presence of a τ had is indicated by the presence of its visible decay products (τ had-vis ).
Events in the τ had τ had channel are selected using a combination of triggers that require either one or two τ had-vis . The triggers used between 2016 and 2018 that require two τ had-vis also require one or two anti-k t R = 0.4 jets. Additionally, events are required to have exactly two τ had-vis with opposite electric charge and exactly two b-tagged anti-k t R = 0.4 jets. Events are rejected if they contain an electron or muon. The invariant mass of the di-τ  [55], which uses both τ had-vis and E miss T with the assumption that E miss T is exclusively from the neutrinos from the τ lepton decays. Events are required to have m MMC ττ > 60 GeV. Events in the τ lep τ had channel are selected using a combination of single-lepton triggers (SLTs) that require a single lepton and lepton-plus-τ had-vis triggers (LTTs) that require a single lepton and a τ had-vis . The LTTs used in 2017 and 2018 one or two jets. Events are split into the SLT and LTT categories based on the triggers used to select them. If an event passes SLT and LTT requirements, it is removed from the LTT category. Events in both categories are required to have exactly one electron or muon, a τ had-vis with opposite electric charge to the lepton, and exactly two b-tagged anti-k t R = 0.4 jets. The MMC is used to estimate the di-τ mass from the τ had-vis , the lepton, and E miss T . Events are required to have m MMC ττ > 60 GeV. Additionally, events are required to have an invariant mass of the two b-tagged jets of m bb < 150 GeV.
The backgrounds in both channels are estimated using a combination of simulation and data-driven methods. The primary backgrounds come from top quark, Z+jets, W+jets, diboson, single Higgs boson, and QCD multijet processes. These backgrounds contain τ had-vis originating from a τ had decay (true τ had-vis ), from a misidentified quark-or gluoninitiated jet (fake τ had-vis ), or from a misidentified lepton. Background events with true τ had-vis and misidentified leptons as well as background with fake τ had-vis other than tt and QCD multijet are estimated using simulation. The normalizations of the tt and Z plus heavy flavor jets (Z+HF) backgrounds are derived using a dedicated control region.
The QCD multijet and tt backgrounds containing one or more fake τ had-vis are estimated using data-driven methods. In the τ lep τ had channel, the QCD multijet and tt backgrounds are estimated using a fake factor method. Fake factors are found separately for each background process using dedicated control regions. In the τ had τ had channel, two different data-driven methods are used for the QCD multijet and tt backgrounds. The QCD multijet background is estimated using a fake factor method similar to the one used in the τ lep τ had channel. The tt background with fake τ had-vis in the τ had τ had channel is estimated using simulation with corrections derived from data in the tt control region from the τ lep τ had data-driven background estimate. The m HH distributions for the τ had τ had and τ lep τ had SLT categories are shown in Figure 19 for the background estimate, data, and several signal hypotheses. The QCD multijet and tt backgrounds containing one or more fake τ had-vis are 522 estimated using data-driven methods. In the τ lep τ had channel, the QCD multijet and tt 523 backgrounds are estimated using a fake factor method. Fake factors are found separately 524 for each background process using dedicated control regions. In the τ had τ had channel, 525 two different data-driven methods are used for the QCD multijet and tt backgrounds.

526
The QCD multijet background is estimated using a fake factor method similar to the 527 one used in the τ lep τ had channel. The tt background with fake τ had-vis in the τ had τ had 528 channel is estimated using simulation with corrections derived from data in the tt control 529 region from the τ lep τ had data-driven background estimate. The m HH distributions for 530 the τ had τ had and τ lep τ had SLT categories are showin in Figure 19 for the background 531 estimate, data and several signal hypotheses.   Figure 19. Distributions of m HH in the (a) τ had τ had and (b) τ lep τ had SLT categories after the background-only fit. The background distribution is shown with its total post-fit uncertainty. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. Select signal hypotheses are overlaid, normalized to a cross-section of 1 pb, as well as the SM HH production hypothesis scaled by a factor of 100. The bottom panel shows the relative difference between the data and background model. [54] Parameterized neural networks (PNNs) are used to distinguish signal from back-533 ground. The PNNs are parameterized in the mass of the resonance hypotheses, providing 534 good sensitivity and continuity over the range of resonance mass hypotheses. The PNNs 535 are trained using the simulated signal as well as the sum of all background estimates 536 normalized to their SM cross-sections. All backgrounds, including those containing 537 fake τ had-vis , are modelled for the training using simulation, except for the QCD multijet 538 background in the τ had τ had channel, which is estimated using the data-driven technique. 539 Figure 19. Distributions of m HH for data, the background estimates, and example signal hypotheses in the (a) τ had τ had and (b) τ lep τ had SLT categories after the background-only fit. The background distribution is shown with its total post-fit uncertainty. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The spin-0 signal hypotheses are normalized to a cross-section of 1 pb, and the SM HH signal hypothesis scaled by a factor of 100. The bottom panels show the ratio of the data to the background estimates [54].

(b)
Parameterized neural networks (PNNs) are used to distinguish signal from background. The PNNs are parameterized in the mass of the resonance hypotheses, providing good sensitivity and continuity over the range of resonance mass hypotheses. The PNNs are trained using the simulated signal as well as the sum of all background estimates normalized to their SM cross-sections. All backgrounds, including those containing fake τ had-vis , are modeled for the training using simulation, except for the QCD multijet background in the τ had τ had channel, which is estimated using the data-driven technique described above. The PNNs are trained using a variety of kinematic variables including m HH , m MMC ττ , m bb , and angular information about all of the decay products. The distributions of PNN outputs in the τ had τ had channel for two signal mass hypotheses are shown in Figure 20.   Figure 20. Distributions of the PNN output scores for (a) m X = 500 GeV in the τ had τ had channel and (b) m X = 1 TeV in the τ had τ had channel, the background estimate and the signal hypothesis. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The signal hypotheses are scaled to the combined expected limit. The lower panels show the ratio between data and the post-fit background and the hatched band shows the total uncertainty on the background estimate. [54] Upper limits on the production cross-section times branching ratio are set at 95% 554 CL using the asymptotic CL s method and are shown in Figure 21. Limits are shown for 555 the τ had τ had and the τ lep τ had channels separately as well as combined. A broad excess 556 is observed in both channels in the signal mass range 700 GeV < m X < 1.2 TeV. The  Obs.
Obs.  The dominant uncertainty is the statistical uncertainty of data. Systematic uncertainties are evaluated on the modeling of the detector response, modeling of the tt, Z+jets and single Higgs boson background simulations, the data driven background estimates, and the choice of PDF and parton shower simulation used for the signal hypotheses.
A simultaneous binned maximum-likelihood fit is performed to the PNN output distributions in the τ had τ had channel and both τ lep τ had categories and to the m distribution in the control region used for the tt and Z+HF normalizations. The tt and Z+HF normalizations are free parameters in the fit. The binnings of the PNN output distributions are chosen to minimize the number of bins while maximizing the expected sensitivity and ensuring the stability of the fit.
The upper limits on the production cross-section times branching ratio are set at 95% CL using the asymptotic CL s method and are shown in Figure 21. Limits are shown for the τ had τ had and the τ lep τ had channels separately as well as combined. A broad excess is observed in both channels in the signal mass range 700 GeV < m X < 1.2 TeV. The most significant combined excess is at a signal hypothesis mass of 1 TeV with a local significance of 3.0 σ and a global significance of 2.0 +0.4 −0.2 σ. A deficit is observed in the τ lep τ had channel in the low-mass region with the largest deficit at a mass of 280 GeV with a local significance of 2.4 σ. Neither the excess at 1 TeV nor the low-mass deficit are reflected in the m HH distributions shown in Figure 19.   Figure 21. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant ggF HH production in the resolved bbτ + τ − channel using the full Run 2 ATLAS dataset. Limits are shown for the τ had τ had and τ lep τ had channels as well as the statistical combination of the channels. The ±1σ and ±2σ uncertainties on the combined expected limits are shown with the colored bands [54].

ATLAS Boosted bbτ + τ − Search
The second analysis [56] is a search for high-mass ggF resonant HH production using a novel reconstruction technique for boosted final states in the τ had τ had channel. The search is interpreted in the context of a generic narrow spin-0 benchmark model. The signal samples are generated with a width of that is much smaller than the detector resolution. Limits are set on the benchmark model in the mass range from 1 to 3 TeV.
A novel technique is used to reconstruct and identify boosted H → τ had τ had decays such that the hadronic decay products overlap in the detector and cannot be distinguished using standard τ had reconstruction. In this method, a di-τ object is reconstructed from an ungroomed anti-k t R = 1.0 jet with p T > 300 GeV. The constituents of the large-R jet are reclustered into anti-k t R = 0.2 subjets, the two p T -leading of which are used to reconstruct the di-τ system. Tracks are matched to the subjets and are referred to as "τ tracks". Other tracks associated to the large-R jet are referred to as "isotracks". All tracks are required to pass criteria on their association to the primary vertex. The charge of each subjet is defined as the sum of the charge of the τ-tracks, and the charge product of the di-τ object is the product of the charges of the two subjets. An illustration of a di-τ object is shown in Figure 22. Requirements are placed on di-τ objects to identify boosted H → τ + τ − decays. Both of the p T -leading subjets are required to have either one or three associated tracks within ∆R < 0.1 to separate τ had decays from hadronic jets. Additionally, a BDT is trained to differentiate between H → τ + τ − and quark-or gluon-initiated jets. The BDT uses information about the τ tracks, isotracks, p T -leading subjets, and large-R jet. A cut is implemented on the BDT output to have an approximately 60% signal efficiency.
Events are selected using triggers that require one high-p T large-R jet. Additionally, events are required to contain at least two large-R jets, at least one of which is a di-τ object. If multiple di-τ objects exist, the one with the highest p T is selected as the H → τ + τ − candidate. The H → τ + τ − candidate is required to have no more than three subjets, and the p T -leading subjets are required to be within ∆R < 0.8 of each other. The two p T -leading subjets are also required to have charges whose product is ±1. Events are also required to have a large-R jet with a mass larger than 50 GeV and separated from the H → τ + τ − candidate by ∆R > 1.0 to reconstruct the boosted H → bb decay. Finally, events have a minimum E miss T requirement and are vetoed if they contain an electron, a muon, or more than one b-tagged large-R jet.
Three signal regions, two control regions, and three validation regions are defined based on the charge product of the H → τ + τ − candidate, the number of b-tagged track jets in the H → bb jet, the |∆φ| between the E miss T and the H → τ + τ − candidate, the mass of the H → bb jet, and the reconstructed mass of the visible HH system m vis HH . The signal regions are defined with a charge product of −1, exactly two b-tagged track jets in the H → bb jet, |∆φ| < 1.0 between E miss T and the H → τ + τ − candidate, and an H → bb jet mass between 60 and 160 GeV. The signal regions are split by m vis HH . For signal hypotheses with m X ≤ 1.6 TeV, no requirement is placed on m vis HH . For signal hypotheses with m X ≥ 1.6 TeV, m vis HH is required to be at least 900 GeV, and for m X ≥ 2.5 TeV, m vis HH is required to be at least 1.2 TeV. The signal hypotheses with m X = 1.6 and 2.5 TeV are each evaluated with both m vis HH requirements that are relevant. Following the signal region requirements, the primary backgrounds come from QCD multijet with hadronic jets misidentified as di-τ objects as well as ZH and Z+jets events with Z → τ + τ − decays. The QCD multijet background is estimated with a data-driven technique that uses a control region in which the di-τ candidate has a charge product of +1 and is validated in two validation regions. The Z+jets background is estimated using simulation with a normalization derived from data in a control region with zero b-tagged track jets and validated in a validation region. An additional normalization is derived for Z+jets events with heavy-flavor jets from data in a control enriched with Z → + − events where is either an electron or a muon. Other background processes including W+jets, diboson, and events containing top quarks are estimated using simulation.
The distribution of m vis HH in the signal region without any requirement on m vis HH is shown in Figure 23 for data, the background estimate, and two signal hypotheses.
The primary systematic uncertainties come from the data-driven background estimates. Additional systematic uncertainties are evaluated on the signal modeling including the choice of PDF and parton shower modeling. Other systematic uncertainties are found to be negligible compared to the statistical uncertainty.
A single bin profile-likelihood fit is conducted for each signal hypothesis. No significant excesses are observed over the predicted background. The upper limits on the production cross-section times branching ratio are set at 95% CL using the CL s method and are shown in Figure 24.

CMS bbτ + τ − Searches
The CMS Collaboration has conducted two searches for ggF resonant HH production in the bbτ + τ − final state using a partial Run 2 dataset of 35.9 fb −1 of data collected in 2016. The first analysis uses resolved topologies to search for resonance masses between 250 and 900 GeV. The second analysis uses boosted topologies to search for resonance masses between 1 and 4 TeV. 800 Figure 23. Distribution of m vis HH in the signal region without the requirement on m vis HH for data, the background estimate, and two example signal hypotheses. The error bar on the data point represents the Poisson uncertainties corresponding to its event yield. The signal hypotheses are normalized to their expected dross-section exclusion limits. The first and last bins contain the underflow and overflow bin entries, respectively. The hatched band shows the statistical plus systematic uncertainty [56].  Figure 24. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant ggF HH production in the boosted bbτ + τ − channel using the full Run 2 ATLAS dataset. The ±1σ and ±2σ uncertainties on the expected limits are shown with the colored bands. The discontinuities at m X = 1.6 and 2.5 TeV come from the different requirements on m vis HH . The signal hypotheses at both points are each evaluated with the m vis HH requirement used on both sides [56].
The first analysis [57] is a search for ggF resonant HH production with a resolved topology in the τ had τ had and τ lep τ had channels. The search is interpreted in the context of a generic narrow spin-0 benchmark model S. The signal samples are generated with a width that is much smaller than the detector resolution. Limits are set on the benchmark model in the mass range from 250 to 900 GeV.
The analysis is split into three channels based on the decay mode of the τ + τ − system. The channel in which both τ leptons decay hadronically (τ had ) is referred to as τ had τ had ; the channels in which one τ lepton decays hadronically and the other decays leptonically are referred to as τ e τ had and τ µ τ had , depending on the flavor of the lepton, or τ lep τ had when discussed inclusively. The presence of a τ had is indicated by the presence of its visible decay products (τ had-vis ).
Events in the τ had τ had channel are selected using a trigger that requires two τ had-vis . Triggers that require a single lepton are used to select events in the τ lep τ had channels. Events in the τ lep τ had channels are required to have one isolated lepton and one τ had-vis , while events in the τ had τ had channel are required to have two τ had-vis . The τ had-vis in all channels are required to pass isolation criteria. Additionally, the visible decay products of the two τ leptons (τ had-vis or lepton) are required to have opposite sign electric charges. Finally, events are required to contain either two anti-k t R = 0.4 jets, at least one of which is b-tagged, or one b-tagged anti-k t R = 0.8 jet with at least two subjets and a mass larger than 30 GeV. If an event has a large-R jet meeting the criteria, it is categorized as boosted; otherwise, it is categorized as resolved. Events in the resolved category are divided into those with exactly one b-tagged jet and those with at least two b-tagged jets and are further split into low-mass and high-mass categories optimized for signal hypotheses with masses less than and greater than 350 GeV.
The invariant masses of the bb system (m bb ) and the τ + τ − system (m ττ ) are reconstructed. In the resolved category, m bb is reconstructed from the two selected jets, and in the boosted category, it is reconstructed from the selected large-R jet. In all channels, m ττ is reconstructed from the τ lepton visible decay products using the SVFit method [58]. Events in the resolved category are required to satisfy requirements that m bb and m ττ are consistent with m H .
The primary background processes are tt, Z/γ * → τ + τ − , and QCD multijet production. The tt, single Higgs boson, W+jets, and diboson backgrounds are estimated using simulation. The Z/γ * → τ + τ − background is estimated using simulation with corrections derived from data in control regions enriched with Z → µ + µ − decays. The QCD multijet background is estimated using data in control regions in which τ lepton visible decay products fail the isolation criteria and events have same-sign pairs of τ had-vis and leptons.
Two BDTs are used to discriminate signal from background in the τ lep τ had channels in the low-and high-mass categories. The BDTs are trained separately for each category using eight variables that consist of angular separations between particles, E miss T , and multiparticle systems as well as the transverse mass of the E miss T plus the τ had-vis or lepton. A kinematic fit [59] is used to reconstruct an invariant mass m KinFit HH by minimizing a function of the b-jet energy, and the p T balance of the selected jets, τ lepton visible decay products, and E miss T . Figure 25 shows the distributions of m KinFit HH for data, background, and signal hypotheses in the low-mass resolved and boosted categories. The dominant systematic uncertainties come from the QCD multijet background estimate, the cross-section of the simulated backgrounds, and the modeling of the detector response for τ leptons. Systematic uncertainties are additionally evaluated on the detector response modeling for other particles, the lepton trigger modeling, and theoretical uncertainties on the signal simulation.
A simultaneous binned maximum-likelihood fit is performed to the m KinFit HH distributions in all of the analysis categories. No significant excesses are observed over the predicted background. Upper limits on the production cross-section times branching ratio are set at 95% CL using the asymptotic CL s method and are shown in Figure 26. The theoretical prediction for a bulk radion with Λ R = 3 TeV and κl = 35 is also shown. The limits are also interpreted as constraints in the m A − tan β plane of the hMSSM as shown in Figure 27.

CMS Boosted bbτ
The second analysis [60] is a search for ggF resonant HH production with a boosted topology using a partial Run 2 dataset of 35.9 fb −1 of data collected in 2016.
The search is interpreted in the context of narrow spin-0 and spin-2 benchmark models. A bulk radion with Λ R = 1 TeV is used as the spin-0 benchmark model, and a Kaluza-Klein graviton withk = κ/M Pl = 0.5 is used as the spin-2 benchmark model. Both benchmark models are generated with a width that is smaller than the experimental resolution. Limits are set on both benchmark models in the mass range from 900 GeV to 4 TeV. The analysis is split into the τ had τ had and τ lep τ had channels, and events are reconstructed using boosted techniques that make use of large-R jets.
A dedicated algorithm is used to reconstruct boosted H → τ + τ − candidates in the τ had τ had and τ lep τ had channels. First, Cambridge-Aachen R = 0.8 jets are used as seeds. The clustering sequence is undone iteratively until two subjets are found that each have a mass less than 2/3 that of the full jet. A τ lepton reconstruction algorithm is used on both subjets, and a BDT is used to identify τ had . If no τ had is identified, the procedure is repeated with anti-k t R = 0.4 jets as seeds. The H → τ + τ − candidate consists of either two τ had or one τ had and a lepton.
Events are selected with a combination of triggers that require large E miss T along with other criteria such as high-p T jets. Additionally, events are required to contain one H → τ + τ − candidate and a H → bb candidate built from an anti-k t R = 0.8 jet with a mass of 105-135 GeV and either one or two b-tagged subjets. Events are also required to have a large E miss T , a di-τ mass estimated with the SVFit algorithm between 50 and 150 GeV, and an invariant mass m X calculated from the H → bb and H → τ + τ − candidates greater than 750 GeV. Finally, events are rejected if they have any b-tagged anti-k t R = 0.4 jets that do not overlap with the H → bb candidate or if the two τ leptons have a ∆R separation greater than 1.5. Following the selection, events are divided into categories based on the number of τ had and the number of b-tagged subjets in the H → bb candidate.
The primary background processes are tt and V+jets. The tt background is estimated using the simulation with a normalization derived from data in a control region with at least one b-tagged anti-k t R = 0.4 jet that does not overlap with the H → bb candidates. The V+jets background is estimated using data in regions in which the H → bb candidate mass is in the range of 30-65 GeV or greater than 135 GeV. The distributions of both background processes are fit with analytic functions that are used to estimate the background contribution in the signal region. A transfer function is used to extrapolate the V+jets estimate to the signal region. The distributions of m X for data, the background estimate, and an example signal hypothesis in the 1 b-tag categories are shown in Figure 28.  Figure 28. Distributions of m X for data, the background estimate and the spin-0 signal hypothesis with m X = 2 TeV and Λ R = 1 TeV in the (a) τ had τ had 1 b-tag and (b) τ lep τ had 1 b-tag categories. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The lower panel shows the difference between data and the background estimate divided by the statistical uncertainty in data. [60] Systematic uncertainties are evaluated on the background estimate, the choice of  Systematic uncertainties are evaluated on the background estimate, the choice of PDF and scales for the signal, and modeling of the detector response. The dominant uncertainties come from the tt and V+jets estimate normalization, the choice of signal PDF and scales, and the modeling of τ lepton identification.
A combined profile-likelihood fit of signal and background is performed in the m X distribution. No significant excess over the background prediction is observed. The upper limits on the production cross-section times branching ratio are set on both models at 95% CL using the CL s method. The limits are shown in Figure 29. The theoretical predictions for a bulk radion with Λ R = 1 TeV and for a Kaluza-Klein graviton withk = κ/M Pl = 0.5 are also shown. The bulk radion model is excluded for masses below 2.7 TeV.

Comparison of bbτ + τ − Searches
Both ATLAS and CMS performed searches with the bbτ + τ − final state using resolved and boosted topologies. The resolved analyses provide excellent sensitivity to resonance masses below 1 TeV. Both resolved analyses make use of the τ lep τ had and τ had τ had channels. Multivariate techniques are used to discriminate signal from background, and data-driven techniques are used to estimate the background processes that include hadronic jets misidentified as τ leptons. A major difference between the ATLAS and CMS analyses is that the CMS analysis uses m HH as the fit variable while the ATLAS analysis uses the PNN output for the fit. The PNN score fit may be the cause of the broad nature of the excess observed around 1 TeV due to the fact that the resonance mass information is not directly used.
The CMS boosted search makes use of the τ lep τ had and τ had τ had channels using large-R jets to reconstruct the H → τ + τ − decay. This results in sensitivity that is competitive at high resonance masses. The ATLAS boosted search uses the τ had τ had channel with dedicated di-τ objects to reconstruct the H → τ + τ − decay. This provides a significantly better selection efficiency at high resonance masses than resolved reconstruction methods; however, the sensitivity is not competitive with other channels. Future iterations of the analysis would benefit from including the τ lep τ had channel as well as optimizing the di-τ identification criteria.

Searches for H H Production in the bbγγ Channel
Searches in the bbγγ channel benefit from the large H → bb branching ratio and the clean signature and excellent mass resolution from the H → γγ decay. The channel has a very low background and provides excellent sensitivity to low-mass resonance hypotheses, below approximately 350 GeV.

CMS bbγγ Search
The CMS Collaboration conducted a search for ggF resonant HH production in the bbγγ channel [61] using a partial Run 2 dataset of 35.9 fb −1 of data collected in 2016.
The search is interpreted in the context of narrow spin-0 and spin-2 benchmark models. Two bulk radion models with Λ R = 2 TeV and Λ R = 3 TeV are used as the spin-0 benchmark models. Two Kaluza-Klein bulk graviton models with κ/M Pl set to 0.5 and 1.0 are used as the spin-2 benchmark models. All benchmark models are generated with a width of 1 MeV. Limits are set on all benchmark models in the mass range from 260 to 900 GeV.
Events are selected using triggers that require two photons that pass isolation requirements and have an invariant mass of m γγ > 90 GeV. The two p T -leading of which are selected to form the H → γγ candidate and are required to pass identification requirements and have an invariant mass of 100 GeV < m γγ < 180 GeV. The two photons are required to have transverse momenta that are at least 33% and 25% of m γγ . Events are also required to have at least two anti-k t R = 0.4 jets, each of which is separated from the nearest selected photon by ∆R > 0.4. If there are more than two jets, the two with the highest b-tagging scores are selected to form the H → bb candidate. The selected jets are required to have an invariant mass of 70 GeV < m jj < 190 GeV.
Two methods are implemented to improve the mass resolution of the reconstructed objects. A multivariate regression technique based on heavy jet fragmentation features and E miss T is used to correct the energy of the selected jets and improve the reconstructed H → bb mass resolution. Additionally, a modified invariant mass of the γγjj system,M X , is used in place of the invariant mass m γγjj . The modified invariant mass is defined as As shown in Figure 30,M X gives an improved mass resolution for resonant signals. The analysis is split into low-and high-mass categories for signal hypothesis with m X ≤ 600 GeV and m X ≥ 600 GeV, respectively, with m X = 600 GeV being used in both categories. A BDT is trained in each category to separate signal from background. The BDTs use the b-tagging information for the selected jets, transverse momentum balance variables of the γγ and jj systems, and three helicity angles of the γγ, jj, and γγjj systems in the Collins-Soper frame of the four-body system [62] as inputs. The BDTs are trained on the ensemble of signal samples in each category as well as background events taken from data in a control region in which one of the two selected photons fails the identification requirements.
The analysis is performed in unique signal regions for each signal hypothesis, defined using aM X window. Within each of these regions, the BDT output is used to define highand medium-purity categories. The criteria for the purity categories for a signal hypothesis depend on the mass category containing the signal hypothesis.
The signal hypotheses are parameterized using fits of a double-sided Crystal Ball function [63], which consists of a Gaussian core and asymmetric power law tails, to both m γγ and m jj in each signal region. The backgrounds containing a single Higgs boson are modeled using functions fit to the m γγ and m jj distributions and are normalized to their SM predictions. A double-sided Crystal Ball function is fit to the background m γγ distributions. Either a double-sided Crystal Ball function or a Bernstein polynomial is fit to the background m jj distributions, depending on the production mechanism. The γγ+jets background is modeled using first and second order Bernstein polynomials fit to data from the control region that are downsampled to match the total number of data events seen in the signal region.
Systematic uncertainties are evaluated on the single Higgs boson background and signal models. The experimental uncertainties include the modeling of the detector response for photons, jets, and b-tagging. Theoretical uncertainties include uncertainties on α S and the choice of PDF as well as cross-section calculations for the single Higgs boson background processes.
The data are interpreted with an unbinned two-dimensional maximum likelihood fit in the m γγ − m jj plane. No significant excesses are observed over the predicted background. The upper limits on the production cross-section times branching ratio are set on both models at 95% CL using the CL s method. The limits are shown in Figure 31. The theoretical predictions for bulk radions with Λ R = 2 TeV and Λ R = 3 TeV and for bulk gravitons with κ/M Pl = 0.5 and κ/M Pl = 1.0 are also shown.

ATLAS bbγγ Search
The ATLAS Collaboration conducted a search for ggF resonant HH production in the bbγγ final state [64] using the full Run 2 dataset. This corresponds to 136 fb −1 of data collected between 2015 and 2018.
The search is interpreted in the context of a generic narrow spin-0 benchmark model. The signal samples are generated with a width of 10 MeV. Limits are set on the benchmark model in the mass range from 251 GeV to 1 TeV.
Events are selected using triggers that require two photons. Additionally, events are required to contain at least two photons that pass isolation requirements, the two p T -leading of which are selected to form the H → γγ candidate. The selected photons are required to have an invariant mass of 105 GeV < m γγ < 160 GeV, and the two photons are required to have transverse momenta that are at least 35% and 25% of m γγ . Events are also required to have exactly two b-tagged anti-k t R = 0.4 jets, fewer than six jets with |η| < 2.5, and exactly zero electrons or muons.
The signal region is defined using a mass window of 120 GeV < m γγ < 130 GeV, which is determined by the H → γγ mass resolution. The sideband region is defined as 105 GeV < m γγ < 160 GeV, excluding the signal region.
A modified invariant mass of the bbγγ system, m * bbγγ , is used in place of the invariant mass m bbγγ . The modified mass is defined as As shown in Figure 32, m * bbγγ gives an improved resolution for resonant signals. A BDT-based method is used to separate signal from background. All signal samples are used together in the training, reweighted event-by-event to have the same m * bbγγ distribution as the background events. Two separate BDTs are trained: one optimized for separating signal from the continuum γγ background and the other optimized for separating signal from backgrounds containing a single Higgs boson (primarily ZH and ttH). The outputs of the two BDTs are combined to construct a single discriminating variable. A minimum threshold on the combined BDT score is chosen for each signal hypothesis to maximize signal significance while ensuring a minimum number of events in the sideband region.
Following the selection requirements, the background consists of the production of single Higgs bosons in the H → γγ decay mode, SM production of HH in the bbγγ decay mode, and continuum γγ+jets events that either contain at least two real photons or one real photon and one or more jets misidentified as photons. The continuum γγ+jets background is parameterized with a functional form fit to a high-statistics simulated background template and normalized using data in the sideband region. The potential bias of various analytic functions on the signal event yield, referred to as "spurious signal", is used to determine the optimal functional form. The functional form with the minimal number of degrees of freedom while remaining below a spurious signal threshold is chosen as exp(a · m γγ ).
The signal as well as backgrounds involving a single Higgs boson and SM HH production are parameterized by fitting the m γγ distribution in simulated samples. The distributions are fit with a double-sided Crystal Ball function [65]. The m γγ distributions for signal and background simulation as well as data after the optimized BDT score requirements are applied are shown for two signal hypotheses in Figure 33. Observed limit (95% CL) Expected limit (95% CL) σ 1 ± Expected limit σ 2 ± Expected limit Figure 34. Expected and observed 95% confidence level limits on the production cross-section The primary systematic uncertainty comes from the spurious signal effects used to select the functional form for the continuum γγ background estimate. Other systematic uncertainties are evaluated on the modeling of the diphoton triggers and the photon and jet modeling, as well as on the PDF and parton showering used in the signal simulation.
A maximum likelihood fit of the m γγ distribution is used to interpret the data. No significant excesses are observed over the predicted background. The upper limits on the production cross-section times branching ratio are set at 95% CL using the CL s method and are shown in Figure 34.

Discussion of bbγγ Searches
ATLAS and CMS both performed searches in the bbγγ channel. The H → γγ decay provides excellent m H resolution and the addition of two b-tagged jets results in a small background.
Both analyses make use of corrected four-body invariant masses to improve the signal mass resolution otherwise degraded by the b-jet energy resolution. Both analyses also use multivariate techniques to discriminate signals from background. Finally, both analyses fit data in the m γγ sidebands to estimate the continuum γγ background in the signal mass window. The ATLAS search fits the m γγ spectrum, while the CMS search does a simultaneous fit to the m γγ and m jj spectra. It is unclear whether the additional constraint of m jj improves the sensitivity due to the poor m jj resolution. Both analyses are dominated by statistical uncertainties but provide excellent sensitivity for low resonance masses. Observed limit (95% CL) Expected limit (95% CL) σ 1 ± Expected limit σ 2 ± Expected limit Figure 34. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant ggF HH production in the bbγγ channel using the full Run 2 ATLAS dataset. The ±1σ and ±2σ uncertainties on the expected limits are shown with the colored bands [64].

Searches for H H Production in the bbV V * Channels
The bbVV * channels offer several advantages to searches for resonant HH production. The large H → bb and H → WW * branching ratios make bbWW * the channel with the second highest HH branching ratio. For both bbWW * and bbZZ * , the final states including one or more leptons helps to significantly reduce QCD multijet backgrounds. Additionally, the presence of intermediate W or Z bosons resonances provides a way to further reduce backgrounds. A complication of these channels is that one of the vector bosons is produced off its mass shell. The major backgrounds for these channels are tt and W+jets.

CMS bbVV * → bb νν Search
The CMS Collaboration has conducted a search for ggF resonant HH production in the bbWW * and bbZZ * channels with the bb νν final state [66] using a partial Run 2 dataset of 35.9 fb −1 of data collected in 2016.
The search is interpreted in the context of narrow spin-0 and spin-2 benchmark models. The RS model in warped extra dimensions with Λ R = 1 TeV and κl = 35 is used for both scenarios. A bulk radion is used as the spin-0 benchmark, and a Kaluza-Klein graviton withk = κ/M Pl = 0.1 is used as the spin-2 benchmark.
Events are selected using triggers that require two leptons. Additionally, events are required to contain two b-tagged anti-k t R = 0.4 jets and two leptons with opposite electric charge and a dilepton invariant mass between 12 GeV and m Z − 15 GeV, where m Z is the mass of the Z boson. The dilepton invariant mass requirement suppresses backgrounds containing quarkonia resonances and Z+jets processes while preserving signal events in the bbWW * channel and in the bbZZ * channel when the off-shell Z boson decays leptonically.
The dominant background is tt production. The second largest background process is Drell-Yan (DY) production. All backgrounds except DY e + e − and µ + µ − production are estimated using simulation. The e + e − and µ + µ − DY backgrounds are estimated using a data-driven method. This is performed using data from a control region that has no b-tagged jet requirement. Events from this region are weighted to estimate the effects of the b-tagging requirement in the signal regions.
A PNN is used to discriminate signal from the large tt background. The PNN makes use of kinematic variables and the flavor of the dilepton pair (e + e − , µ + µ − , or e ± µ ∓ ). Three regions are defined using the invariant mass of the two b-tagged jets, m jj . A signal-enriched region is defined as m jj ∈ [75, 140) GeV, and background-enriched regions are defined on either side of this range. The PNN output is used as the final discriminant in each region.
Systematic uncertainties are evaluated on the data-driven background estimate, theoretical predictions, and detector response. Theoretical uncertainties are evaluated on the choice of PDF as well as the renormalization and factorization scales used in simulation. Detector response uncertainties are evaluated on trigger efficiencies, lepton identification, jet energy measurements, and b-tagging. The data-driven background uncertainty consists of propagated detector response uncertainties plus an additional normalization uncertainty.
Results are obtained through a maximum likelihood fit of the PNN discriminant in the three m jj regions in each of the three lepton flavor categories. No significant excesses are observed over the predicted background. The upper limits are set at 95% CL using the asymptotic CL s method on the production cross-section times branching ratio for both benchmark models.The upper limits on cross-section times branching ratio are shown in Figure 35. The theoretical predictions for a bulk radion with Λ R = 3 TeV and κl = 35, for a bulk gravitons with κ/M Pl = 0.1 and κl = 35 are also shown.

CMS bbZZ * Search
The CMS Collaboration conducted a search for ggF resonant HH production in the bbZZ * channel with the bb final state [67] using a partial Run 2 dataset of 35.9 fb −1 of data collected in 2016.
The search is interpreted in the context of narrow spin-0 and spin-2 benchmark models. The RS model in warped extra dimensions with Λ R = 1 TeV and κl = 35 is used for both scenarios. A bulk radion is used as a spin-0 benchmark, and a Kaluza-Klein graviton with k = κ/M Pl = 0.1 is used as a spin-2 benchmark. In addition to setting cross-section times branching ratio limits on these, scans of the Λ R and κ/M Pl parameters are performed. The search is also interpreted in the context of the N2HDM where the masses of H 1 and H 2 are both set to 125 GeV and H 3 is treated as the heavy spin-0 resonance. Two benchmark points are chosen with tan β = 0.5 and 2.0. In both cases, the scalar vev is set to 45 GeV, and the mixing angles α 1 , α 2 and α 3 are set to 076, 0.48, and 1.00, respectively. This results in an increase of the HH branching fraction to bbZZ * of 33% relative to the SM prediction for tan β = 0.5 and 5% for tan β = 2.0. All of the signal samples are simulated in the mass range of 260 ≤ m X ≤ 1000 GeV with narrow resonances.
The analysis is split into channels to separate the bb qq (denoted as the bb jj channel) and bb νν decay modes. Events in both channels are selected using triggers that require two muons or two electrons. Events in both channels are also required to have two sameflavor leptons with opposite electric charge. The hadronic decay products in both channels are reconstructed using anti-k t R = 0.4 jets. A BDT is trained for each mass and spin hypothesis in each lepton channel to maximize sensitivity, for a total of 48 BDTs.
Events in the bb νν channel are required to have a dilepton invariant mass of m ≥ 15 GeV and at least two jets, at least one of which is b-tagged. Additionally, events are required to have a minimum E miss T , the threshold of which is dependent on m X . The H → bb candidate is reconstructed from the pair of jets with the highest values of the b-tagging discriminant if there are at least two b-tagged jets. If there is a single b-tagged jet, it is paired with the jet that gives an invariant mass closest to 125 GeV to form the H → bbcandidate. The H → ZZ * candidate is reconstructed from the pair of leptons and the E miss T . A BDT is trained for each signal hypothesis with 22 kinematic variables as inputs and events are required to pass a minimum threshold of the BDT score.
Events in the bb jj channel are required to have m ≥ 15 GeV and at least four jets, at least one of which is b-tagged. Additionally, events are required to fail the m X -dependent E miss T criteria applied in the bb νν channel. The H → bb candidate is reconstructed following the same method as in the bb νν channel. The H → ZZ * candidate is reconstructed from the two leptons and the two remaining jets that result in the jj invariant mass that is closest to 125 GeV. A BDT is trained for each signal hypothesis with nine kinematic variables as inputs.
The main background processes in the bb νν channel are tt and Z/γ * +jets, with smaller contributions from other non-QCD SM processes. The QCD multijet background is negligible. The non-negligible backgrounds are estimated using simulation. The normalization for the tt and Z/γ * +jets backgrounds are derived from data in dedicated control regions that are fit simultaneously with the signal region in the distribution of the pseudo transverse mass of the HH systemM T (HH).
The main background processes in the bb jj channel are tt and Z/γ * +jets, with smaller contributions from other SM processes, including QCD multijet events. The non-QCD backgrounds are all estimated using simulation. The tt and Z/γ * +jets backgrounds are normalized using data in dedicated control regions. The QCD multijet background is estimated using data events with leptons with the same sign electric charge.
Systematic uncertainties are evaluated on the background model, theoretical predictions, and detector response. The background modeling systematic uncertainties include the normalization and shape determinations. The theoretical uncertainties include the normalization and factorization scales of the background simulation as well as the PDF and α S for the signal simulation. The dominant detector systematic uncertainty is related to the energy of jets.
Results are obtained through a binned maximum likelihood fit. In the bb jj channel, this is performed on distributions of the BDT scores. In the bb νν channel, theM T (HH) is fit simultaneously in the signal region and both control regions for each signal hypothesis.
No significant excesses are observed over the predicted background. The upper limits are set at 95% CL using the asymptotic CL s method on the production cross-section times branching ratio, as well as on the Λ R and thek = κ/M Pl parameters of the RS model. Limits are set in each channel separately as well as in both channels combined statistically. The upper limits on cross-section times branching ratio are shown in Figure 36 for both channels as well as the combination. The theoretical predictions for a bulk radion with Λ R = 3 TeV and κl = 35, for a bulk gravitons with κ/M Pl = 0.1, and for N2HDM models with tan β = 0.5 and 2.0 are also shown. Combined exclusion limits on the RS model parameters are shown in Figure 37 along with limits from each channel individually.

CMS bbWW * Search
The CMS Collaboration has conducted a search for ggF resonant HH production in the bbWW * and bbτ + τ − channels in the bb+leptons final state [68] using the full Run 2 dataset of 138 fb −1 of data collected between 2015 and 2018. The search includes the single lepton final state from the HH → bbWW * → bb νqq decay mode and two lepton final states from the HH → bbWW * → bb ν ν and HH → bbτ + τ − → bb νν νν decay modes.
The search is interpreted in the context of narrow spin-0 and spin-2 benchmark models. A bulk radion with Λ R = 3 TeV and a branching fraction to HH of 25% is used as the spin-0 benchmark model. Two Kaluza-Klein bulk graviton models with branching fractions to HH of 10% andk = κ/M Pl set to 0.3 and 0.5 are used as the spin-2 benchmark models. Both benchmark models are generated with a width of 1 MeV. Limits are set on all benchmark models in the mass range from 800 GeV to 4.5 TeV.
The analysis is split into single lepton (SL) and dilepton (DL) channels. These channels are further split into a total of twelve analysis categories, based on the lepton multiplicity, lepton flavor, b-tagging information, and kinematics of the leptonic decay products. The H → bb decay is reconstructed using a single large-R jet in all categories to exploit the boosted topology.
Events are selected with a combination of triggers that require a single electron or muon (denoted as leptons) or a large scalar sum of hadronic jet p T . After the trigger selection, events are required to contain at least one lepton and a large scalar sum of jet p T . Events in the DL channel are required to have exactly two leptons with opposite electric charge that pass isolation criteria optimized for nearby leptons in boosted topologies and pass loose impact parameter requirements to accept leptons from H → τ + τ − decays. Events in the SL channel are required to have exactly one lepton that passes isolation criteria optimized for topolgies without nearby leptons and no more than one lepton that passes the DL lepton criteria. Events in both channels are required to contain exactly one anti-k t R = 0.8 jet, denoted as the bb jet, that is tagged as being consistent with a H → bb decay using the DeepAK8 tagger. Events with a b-tagged anti-k t R = 0.4 jet that is ∆R > 1.2 from the bb jet are vetoed to reduce the tt background.
Events in the SL channel are required to have a second anti-k t R = 0.8 jet, denoted as the qq jet, that is consistent with a W ( * ) → qq decay. The qq jet is chosen as the anti-k t R = 0.8 jet that is closest in ∆R to the lepton. The qq jet is required to be within ∆R < 1.2 of the lepton and to have substructure that is consistent with a two-body decay. The bb jet is required to be ∆R > 2.0 from the lepton and ∆R > 1.6 from the qq jet. In DL channel events, the bb jet is required to be ∆R > 2.0 from the two lepton system ( ) and neither lepton is allowed to overlap with the bb jet cone.
A likelihood-based method is used to fully reconstruct the H → WW * decay in the SL channel. A likelihood function is constructed from five probability density functions to extract the three components of the neutrino momentum, an indicator of whether the hadronically decaying W boson is on-or off-shell, and a correction factor for the p T of the qq jet. The full H → WW * system is reconstructed from the lepton, the extracted neutrino momentum, and the corrected qq jet. Additionally, a discriminating variable D νqq is constructed using the likelihood function.
The H → WW * decay is reconstructed in the DL channel with assumptions about the decay kinematics. The η of the two neutrino system νν is set to that of . Additionally, the νν invariant mass is set to 55 GeV, which is the average of the truth-level distribution in simulated signals. The H → WW * system is reconstructed from νν and .
In both channels, the HH invariant mass m HH is calculated from the reconstructed H → WW * and the bb jet.
Additional criteria are used to separate signal from background. In the SL channel, the ratio of the H → WW * p T to m HH , denoted as p T /m, is required to be greater than 0.3. In the DL channel, several criteria are applied to and E miss T . The system is required to have an invariant mass between 6 and 75 GeV and a angular separation of ∆R < 1.0. Additionally, E miss T and are required to have a ∆φ separation less than π/2. Finally, events are required to have E miss T > 85 GeV. Events are split into twelve categories based on several properties. The lepton flavor and the DeepAK8 discriminant D Z/H→bb on the bb jet are used in both the SL and DL channels. Additionally, discriminants on the H → WW * system are used to create highand low-purity categories in the SL channel. These discriminants are D νqq and the τ 2 /τ 1 N-subjettiness ratio [69] for the qq jet, which is a substructure variable used to identify jets that are consistent with a two-body hadronic decay. The category labels are defined in Tables 1 and 2. Two Control Regions (CRs) are defined to validate the SM background estimation and systematic uncertainties. The top CR is defined to select events containing top quarks, especially tt. This CR is defined by requiring at least one b-tagged anti-k t R = 0.4 jet that is far from the bb jet. Additionally, requirement on p T /m is removed for the SL channel, and the ∆R requirement is changed to ∆R > 0.4 for the DL channel. The top CR is divided into the same twelve categories as the signal region. The nontop CR is defined to select events from Z/γ * +jets, W+jets and QCD multijet processes. This CR is defined in the same way as the signal region, except the D Z/H→bb requirement on the bb jet is changed to reject events containing H → bb decays. The nontop CR is divided into six categories following the lepton flavor and H → WW * purity discriminants. The background is split into components using generator-level information to ensure distinct m bb distributions. These components are defined by the number of generator-level quarks from the immediate decay of a top quark or vector boson that are within ∆R < 0.8 of the reconstructed bb jet axis. The "m t " component consists of events in which three quarks fulfill the criterion. The "m W " component consists of events in which two quarks do. The "lost t/W" component consists of events in which one quark does. The "q/g" component consists of all other events. The first three consist primarily of tt events, and the q/g component consists primarily of W+jets and QCD multijet events.
Signal and background are modeled in the m HH − m bb plane in each of the twelve event categories using smooth templates created from simulation. The background templates are modeled using conditional probabilities of m bb as a function of m HH based on distributions of the two variables in simulation. The signal templates are modeled similarly but using conditional probabilities of m bb and m HH as a function of the resonance mass m X . The m bb and m HH distributions in the µ bT LP category are shown in Figure 38 for the background estimate, data, and example signal models. Alternate background templates for all four components are built with parameters proportional to m HH and 1/m HH . For the "q/g" and "lost t/W" background components, additional templates are built with parameters proportional to m bb and 1/m bb .
Systematic uncertainties are evaluated on the background and signal models. The background uncertainties include variations of the background component normalizations, the bb jet mass, and the template definitions. The signal uncertainties include uncertainties on theoretical predictions, luminosity measurements, pileup modeling, and detector response modeling for leptons, jet mass and energy scale, jet substructure, and b-tagging.
The data are interpreted with a simultaneous two-dimensional maximum likelihood fit in the m bb − m HH plane in all twelve categories. No significant excesses are observed over the predicted background. The upper limits on the production cross-section times branching ratio are set on both models at 95% CL using the CL s method and are shown in Figure 39 H → WW * system is reconstructed from the lepton, the neutrino solution and W had .

1148
Three signal regions are defined to be sensitive to different ranges of resonance 1149 masses. The first targets signal hypotheses with a mass of 500 GeV. The second is the 1150 Figure 38. Distributions of (a) m bb and (b) m HH for data, the fitted background estimate, and example signal models normalized to a cross-section times branching ratio of 0.2. The lower panels show the ratio of data to the background estimate. The hatched band shows the background shape uncertainty from the 2D maximum likelihood fit to data in the m bb − m HH plane. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields [68].

ATLAS bbWW * Search
The ATLAS Collaboration has conducted a search for ggF resonant HH production in the bbWW * channel in the bb νqq decay mode [70] using a partial Run 2 dataset of 36.1 fb −1 of data collected in 2015 and 2016.
The search is interpreted in the context of narrow spin-0 and spin-2 benchmark models. A generic scalar is used as the spin-0 benchmark model. Two Kaluza-Klein bulk graviton models with c = κ/M Pl set to 1.0 and 2.0 are used as the spin-2 benchmark models. The benchmark models are generated with a narrow width that is much smaller than the detector resolution. Limits are set on all benchmark models in the mass range from 500 GeV to 3 TeV.
The analysis is split into resolved and boosted channels that are sensitive to resonance masses less than and greater than approximately 1 TeV, respectively. The resolved analysis uses anti-k t R = 0.4 jets to reconstruct all hadronic decay products. The boosted analysis uses a combination of anti-k t R = 0.4 jets and anti-k t R = 1.0 (large-R) jets.
Events in both channels are selected with a combination of triggers that require a single lepton. Events are also required to contain at least one lepton that fulfills criteria on the impact parameters relative to the primary vertex.
Events in the resolved channel are required to have at least four anti-k t R = 0.4 jets, exactly two of which are required to be b-tagged. The two b-tagged jets are used to reconstruct the H → bb candidate. If there are exactly two untagged jets, they are used to reconstruct the hadronic W boson candidate W had . If there are more than two untagged jets, the three p T -leading jets are considered, and the pair with the smallest ∆R separation is used to reconstruct W had .
The full H → WW * → νqq topology is reconstructed by solving for the z-component of the neutrino momentum using the Higgs boson mass constraint. If two real solutions exist, the result that minimizes ∆R between the neutrino and the lepton is used. If only a complex solution exists, the real part of the solution is used. The H → WW * system is reconstructed from the lepton, the neutrino solution, and W had .
Three signal regions are defined to be sensitive to different ranges of resonance masses. The first targets signal hypotheses with a mass of 500 GeV. The second is the low-mass range that targets resonance masses between 600 GeV and 1.3 TeV. The third is the high-mass range that targets resonance masses between 1.4 and 3 TeV.
A set of requirements on kinematic variables is optimized to maximize the expected sensitivity to signal hypotheses in each signal region. The variables used are E miss T , the mass and p T of the WW * system, the mass and p T of the bb system, and the ∆R separation between the two reconstructed W bosons. The same selection is used for spin-0 and spin-2 signal hypotheses.
The dominant backgrounds come from tt events with subdominant contributions from QCD multijet, W/Z+jets, single top quark, and diboson processes. The tt background is estimated using simulation and then normalized using data in tt-enriched control regions that correspond to each signal region. The QCD multijet background is estimated in each signal and tt control region using a data-driven method. The other backgrounds are all estimated using simulation.
The m HH distributions in the low-mass and high-mass regions are shown in Figure 41 for data, the background estimate, and example signal hypotheses. Figure 40. Expected 95% CL upper limits on the production cross-section times branching ratio of each of the twelve event categories individually. [68] low-mass range that targets resonance masses between 600 GeV and 1.3 TeV. The third is 1151 the high-mass range that targets resonance masses between 1.4 TeV and 3 TeV.  The m HH distributions for the low-mass and high-mass regions are shown in Figure   1164 41 for the background estimate, data, and example signal hypotheses.  Systematic uncertainties are evaluated on the background and signal models as well as detector response modeling. The background uncertainties include components for the tt simulation and normalization method, the W+jets scale, PDF, and generator modeling, as well as the data-driven QCD multijet estimate. Systematic uncertainties are derived for the scale, PDF set, and parton shower used for the signal samples. Detector systematic uncertainties are used for all physics objects and are dominated by uncertainties on the energy scale of jets.
A maximum likelihood fit is performed on the m HH distribution simultaneously in the signal and control regions. The tt and QCD multijet normalizations are allowed to float in the fit. No significant excess are observed over the predicted background. The upper limits on the production cross-section times branching ratio are set on the signal models at 95% CL using the CL s method across the full resonance mass range. The limits are shown for each signal model in Figure 42. Events in the boosted channel are required to have at least one large-R jet with an angular distance of ∆R > 1.0 from the lepton. The large-R jet with the highest p T is used as the H → bb candidate. The H → bb candidate is required to have a mass between 30 and 300 GeV. Events are required to have at least two R = 0.4 jets with an angular distance of ∆R > 1.4 from the H → bb candidate. The hadronically and leptonically decaying W bosons are reconstructed using the same method as the resolved channel. Events are rejected if any jets with ∆R > 1.4 from the H → bb candidate are b-tagged.
The signal region is defined by requiring the H → bb candidate to have at least two associated track jets and that the two p T -leading associated track jets be b-tagged. Additionally, the H → bb candidate is required to have a mass between 90 and 140 GeV. A minimum E miss T requirement is applied to reduce the QCD multijet background. A validation region is defined to assess the modeling of the tt background. This region is defined by requiring the H → bb candidate have a mass that lies outside the range of 90-140 GeV.
As in the resolved channel, the dominant backgrounds come from tt events with subdominant contributions from QCD multijet, W/Z+jets, single top quark, and diboson processes. The non-QCD backgrounds are all estimated using simulation. The QCD multijet background is estimated using a data-driven technique that is very similar to the one used in the resolved channel.
The m HH distribution for data, the background estimate, and example signal hypotheses in the boosted channel after the global likelihood fit is shown in Figure 43.  The systematic uncertainties in the boosted channel are evaluated in a very similar way to those in the resolved channel. Additional detector response uncertainties are evaluated on the large-R jet mass and energy as well as b-tagged track jets.
A maximum likelihood fit is performed on the m HH distribution in the signal region. Unlike in the resolved analysis, the QCD multijet estimate is not allowed to float in the fit. No significant excess are observed over the predicted background. The upper limits on the production cross-section times branching ratio are set on the signal models at 95% CL using the CL s method across the full resonance mass range. The limits are shown for each signal model in Figure 44. The resolved and boosted channels are both used to set limits across the full resonance mass range of 500 GeV to 3 TeV. The final analysis limits are found by using the limits from whichever channel are the strongest at each resonance mass hypothesis as shown in Figure 45. This results in a discontinuity in the limits at 1.3 TeV for the spin-0 benchmark model and at 800 GeV for the spin-2 benchmark models.

Discussion of bbVV * Searches
ATLAS and CMS have both performed searches in bbVV * channels. The large H → bb and H → VV * branching fractions and the good QCD multijet suppression due to leptonic final states help to make these channels sensitive across a wide range of resonance masses.
CMS performed two analyses that search for low-mass resonances in the 2-lepton bbVV * channel. The first is optimized for the bbWW * channel, and the second is optimized for the bbZZ * channel in both the bb qq and bb νν decay modes. Both of these channels have complex signature, and machine learning techniques are used to maximize sensitivity. These channels are nonetheless limited in their sensitivity compared to other channels at low mass.
The ATLAS bbWW * analysis uses resolved techniques to identify and reconstruct the 1-lepton H → WW * system. This is complicated by the fact that half of the signal events contain a very soft lepton from the decay of the off-shell W boson, decreasing the efficiency of the single lepton triggers. The resolved analysis is further complicated by the fact that in the other half of signal events, the hadronic W boson is off-shell, resulting in low-p T jets. Additionally, the neutrino reconstruction method is detrimental to the analysis sensitivity. In a significant fraction of events, using the H mass constraint results in the correct solution being mathematically impossible to obtain. In addition, the H mass constraint strongly sculpts the background to look like signal. In the boosted analysis, the hadronic W boson decay products overlap with each other and with the lepton for the resonance masses that result in a boosted H → bb decay. As a result, the overlap removal between the lepton and jets results in a low signal selection efficiency and inclusion of radiated jets in the H → WW * reconstruction. Future iterations of the analysis would benefit from additional triggers. The resolved analysis would benefit from more sophisticated selection and reconstruction techniques. The boosted analysis would benefit significantly from particle flow large-R jets and leptons for the H → WW * reconstruction.
The CMS boosted bbWW * analysis makes use of particle flow large-R jets and demonstrates the sensitivity of the 1-and 2-lepton channels for large resonance masses. For masses above 3 TeV, the single muon channel is competitive with the bbbb channel. For masses above 1 TeV, the 2-lepton channel can provide additional sensitivity that is comparable to the 1-lepton and bbbb channels.

Searches for H H Production in Other Channels
Other HH decay channels that do not include an H → bb decay offer additional ways to search for signals. The relatively large H → WW * branching fraction is exploited along with the low backgrounds and clean signatures provided by diphoton and multilepton final states.

ATLAS WW * WW * Search
The ATLAS Collaboration has conducted a search for ggF resonant HH production in the WW * WW * channel in the multilepton final state [71] using a partial Run 2 dataset of 36.1 fb −1 of data collected in 2015 and 2016.
The search is interpreted in the context of a generic narrow spin-0 benchmark model. The signal samples are generated with a width that is much smaller than the detector resolution. Limits are set on the benchmark model in the mass range from 260 to 500 GeV.
The analysis is split into channels by the number of leptons in the final state. The channels are defined as having two leptons with the same sign electric charge, three leptons, and four leptons. The two-lepton channel is split into three categories based on the flavor of the leptons (ee, eµ and µµ). The three-lepton channel is split into two categories based on the number of same-flavor opposite-sign (SFOS) lepton pairs (N SFOS = 0 and N SFOS = 1, 2). The four-lepton channel is split into four categories based on the four lepton invariant mass (m 3 < 180 GeV and m 3 > 180 GeV) and the number of SFOS lepton pairs (N SFOS = 0, 1 and N SFOS = 2).
Events are selected using triggers that require either one or two leptons. In the two lepton channel, events are required to have exactly two leptons with the same sign electric charge and at least three anti-k t R = 0.4 jets. Events in the three lepton channel are required to have exactly three leptons with a summed electric charge of ±1 and at least two jets. Both channels have additional requirements on E miss T and the p T of the leptons. Events are rejected if they contain a SFOS lepton pair with an invariant mass m of |m − m Z | < 10 GeV or m < 15 GeV. In the four lepton channel, events are required to have exactly four leptons with a summed electric charge of 0. Events are rejected if they have a SFOS lepton pair with m < 4 GeV. Events in all channels are rejected if they contain any b-tagged jets.
The backgrounds include processes containing prompt leptons, leptons with incorrectly identified charges, and fake leptons that consist of nonprompt leptons and hadronic jets misidentified as leptons. The prompt lepton backgrounds include diboson, triboson, and ttZ processes and are estimated using simulation. The two-and three-lepton channels have nonprompt lepton backgrounds that come from the conversion of prompt photons in Vγ processes that are estimated using simulation. The remaining fake lepton background contributions and the incorrectly identified charge background are estimated using datadriven methods. Figure 46 shows the invariant masses of the two-and four-lepton systems in the corresponding channels. Optimized selection criteria consisting of requirements on four variables are used to distinguish signal from background. The set of variables and the requirements on each are optimized independently for each signal hypothesis in the two-and three-lepton channels and for each category in the four-lepton channel. The variables include invariant masses and angular separations of combinations of leptons and jets. The yields in each category after the selection optimized for the m X = 260 GeV signal hypothesis are shown in Figure 47.
Systematic uncertainties are evaluated on the modeling of leptons, jets, and E miss T as well as the data-driven background estimate and theoretical modeling of all simulated samples. The dominant uncertainties come from the jet modeling and the fake lepton background estimate. A statistical analysis using a profile-likelihood-ratio test statistic based on the yields in each analysis category is performed in each channel separately as well as a combination of the three channels. No significant excesses are observed over the predicted background. The upper limits on the production cross-section times branching ratio are set at 95% CL using the asymptotic CL s method and are shown in Figure 48.

ATLAS WW * γγ Search
The ATLAS Collaboration conducted a search for ggF resonant HH production in the WW * γγ channel in the γγ νjj final state [72] using a partial Run 2 dataset of 36.1 fb −1 of data collected in 2015 and 2016.
The search is interpreted in the context of a generic narrow spin-0 benchmark model. The signal samples are generated with a width that is much smaller than the detector resolution. Limits are set on the benchmark model in the mass range from 260 to 500 GeV.
Events are selected using triggers that require two photons. Additionally, events are required to contain at least two photons, the two E T -leading of which are selected to form the H → γγ candidate. The selected photons are required to have an invariant mass of 105 GeV < m γγ < 160 GeV, and the two photons are required to have transverse energies that are at least 35% and 25% of m γγ . Events are also required to have at least two anti-k t R = 0.4 jets, at least one lepton and no b-tagged anti-k t R = 0.4 jets. Finally, events are required to have a large transverse momentum of the H → γγ candidate p γγ T for signal hypothesis with m X > 400 GeV, defining two signal regions.
The signal regions are both defined using a window in the diphoton invariant mass of 121.7 GeV < m γγ < 128.5 GeV, which is determined by the H → γγ mass resolution. The sideband regions are defined as 105 GeV < m γγ < 160 GeV, excluding the signal regions.
The background consists of the production of single Higgs bosons in the H → γγ decay mode, SM production of HH in the WW * γγ decay mode, and continuum γγ+jets events that either contain at least two real photons or one real photon and one or more jets misidentified as photons. The continuum γγ+jets background is parameterized with a second-order polynomial fit to data in the sideband.
The signal as well as backgrounds involving a single Higgs boson and SM HH production are parameterized by fitting the m γγ distribution in simulated samples. The distributions are fit with a double-sided Crystal Ball function. The m γγ distributions for data and the parameterized background fits are shown in Figure 49.
The signal samples are generated with a width that is much smaller than the detector 1308 resolution. Limits are set on the benchmark model in the mass range from 260 GeV to 1309 500 GeV.

1310
Events are selected using triggers that require two photons. Additionally, events 1311 are required to contain at least two photons, the two E T -leading of which are selected  Figure 49. Distribution of m γγ for data and the parameterized background fits for resonance masses (a) below 400 GeV and (b) above 400 GeV. The error bars on the data points represent the Poisson uncertainties corresponding to their event yields. The fits are found using the signal-plusbackground likelihood fit. [72] The primary systematic uncertainty comes from the spurious signal effects used The primary systematic uncertainty comes from the spurious signal effects used to select the functional form for the continuum γγ background estimate, similar to the method discussed in Section 6.2. Other systematic uncertainties are evaluated on the modeling of the diphoton triggers and the photon, jet, and b-tagging modeling as well as on the PDF and parton showering used in the signal simulation.
A maximum likelihood fit of the m γγ distribution is used to perform the background estimate and interpret the data. No significant excesses are observed over the predicted background. The upper limits on the production cross-section times branching ratio are set at 95% CL using the CL s method and are shown in Figure 50.

Discussion of Searches in Other Channels
ATLAS has performed searches in the WW * WW * and the WW * γγ channels. Both searches make use of the low backgrounds and clean events offered by H → γγ decays and multilepton final states. The relatively large H → WW * branching ratio is offset by the small branching ratios of W boson decays to leptons. Both searches are performed for resonance masses below 500 GeV where they are expected to be most sensitive. The low branching ratios, poor mass resolutions due to neutrinos, and the very low statistics result in these channels not being as sensitive as channels that include H → bb decays. Future datasets that are significantly larger may make these channels more competitive.

Combination of H H Search Channels
No single decay channel is optimally sensitive to resonant HH production, and different ranges of resonance mass hypotheses are dominated by different channels. A statistical combination of the individual channels is performed to set limits on the signal hypotheses.

ATLAS Channel Combinations
The ATLAS Collaboration has performed two ggF resonant HH channel combinations: the first as a combination of individual channel searches using a partial Run 2 dataset of up to 36.1 fb −1 [73] and the second as a combination of individual channels using the full Run 2 dataset of 126-139 fb −1 [74]. A statistical combination is performed on the results of the individual analyses with a simultaneous fit to the data for the signal hypothesis cross-sections. Nuisance parameters that encode the statistical and systematic uncertainties for each individual channel are included in the fit. The upper limits on the production crosssection times branching ratio are set on the signal models at 95% CL using the asymptotic CL s method. The signal regions used in the simultaneous fit are either orthogonal by construction or have negligible overlap.
The combined results are interpreted in the context of spin-0 and spin-2 benchmark models. A generic scalar S is used as the spin-0 benchmark and is generated with a narrow width that is much smaller than the detector resolution. The spin-2 benchmark are Kaluza-Klein gravitons G KK with κ/M Pl = 1 and 2 with generated widths ranging from 3% to 25% of the resonance masses. Limits are set on all three benchmark models for a mass range from 260 GeV to 3 TeV. The limits for all three benchmark models are shown in Figure 51 with G KK theory predictions overlaid on the respective plots. The G KK model with κ/M Pl = 1 is excluded for the mass range of 310-1380 GeV, and the model with κ/M Pl = 2 is excluded for the mass range of 260-1760 GeV. partial dataset bbbb, bbτ + τ − and bbγγ analyses are similar to the full dataset analyses de-1369 scribed in Sections 4.1, 5.1 and 6.2, respectively. The partial dataset bbWW * , WW * WW * 1370 and WW * γγ analyses are described in Sections 7.4, 8.1 and 8.2, respectively. 1371 The combined results are interpreted in the context of spin-0 and spin-2 benchmark 1372 models. A generic scalar S is used as the spin-0 benchmark and is generated with a 1373 narrow width that is much smaller than the detector resolution. The spin-2 benchmark   Exclusion limits are set on the EWK-singlet model and the hMSSM using a combination of the bbbb, bbτ + τ − , and bbγγ channels. Experimental limits on the spin-0 benchmark model are interpreted as constraints in the m S − m α plane in the EWK-singlet model for tan β = 1 and tan β = 2 as shown in Figure 52. Limits are also set in the sin α − tan β plane for m S = 260 GeV as shown in Figure 53a. The limits are additionally interpreted as constraints in the m A − tan β plane of the hMSSM as shown in Figure 53b.  [77]. The dotted lines show the separation between regions where the width of the resonance is larger than 2% and 5% than its mass. The hatched areas show regions that cannot be excluded because the width of the resonance exceed 10%, which is the maximum mass resolution of the channels included in the combination.   [77]. The dotted lines show the separation between regions where the width of the resonance is larger than 2% and 5% of its mass. The hatched areas show regions that cannot be excluded because the width of the resonance exceeds 10%, which is the maximum mass resolution of the channels included in the combination [73].  Figure 52. Excluded regions of the EWK-singlet model for (a) tan β = 1 and (b) tan β = 2 using the combined upper limits from the spin-0 resonance search. The horizontal lines show indirect constraints from SM Higgs coupling measurements [77]. The dotted lines show the separation between regions where the width of the resonance is larger than 2% and 5% than its mass. The hatched areas show regions that cannot be excluded because the width of the resonance exceed 10%, which is the maximum mass resolution of the channels included in the combination. The full Run 2 dataset combination uses searches in the bbbb [42], bbτ + τ − [54], and bbγγ [64] analyses described in Sections 4.1, 5.1, and 6.2, respectively.
The combined results are interpreted in the context of a spin-0 benchmark model. A generic scalar X is used as the benchmark and is generated with a width of 10 MeV. Limits are set for a mass range from 251 GeV to 3 TeV and are shown in Figure 54. The bbγγ analysis is the most sensitive in the 251-400 GeV mass range, the bbτ + τ − analysis is the most sensitive for 400-800 GeV, and the bbbb analysis is the most sensitive for 800-3000 GeV.

HH) [fb]
ATLAS Preliminary s = 13 TeV, 126 139 fb 1 Spin-0 Observed limit (95% CL) Expected limit (95% CL) Comb. exp. limit ± 1 Comb. exp. limit ± 2 bbbb bb + bb Combined Figure 54. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant ggF HH production combining channels using the full Run 2 ATLAS dataset in the context of a generic scalar X. Limits are shown for the individual channels as well as the statistical combination of the channels. The ±1σ and ±2σ uncertainties on the combined expected limits are shown with the colored bands [74].
The largest deviation from the SM expectation is observed at m X = 1.1 TeV, corresponding to the deviation observed in the bbτ + τ − channel. The local significance of this deviation is 3.2 σ and the global significance is 2.1 σ as shown in Figure 55.  Figure 56 shows an overlay of the limits from the partial Run 2 combination and the limits from the full Run 2 analyses in the bbγγ (Section 6.2), resolved and boosted bbτ + τ − (Section 5.1), and bbbb (Section 4.1) channels.

CMS Channel Combinations
The CMS Collaboration performed a ggF resonant HH channel combination using a partial Run 2 dataset of 35.9 fb −1 [79].
Likelihood fits are performed to combine the results of the individual analyses. Systematic uncertainties from each channel are accounted for as nuisance parameters in the fits, some of which are correlated across multiple channels. The upper limits on the production cross-section times branching ratio are set on narrow-width spin-0 and spin-2 benchmark signal models at 95% CL using the CL s method in the mass range 250 GeV to 3 TeV. The combination is conducted with searches performed in the bbbb [48,80,81], bbτ + τ − [57], bbγγ [61], and bbVV * [66] channels.
The partial dataset boosted bbbb analysis is similar to the full dataset analysis described in Section 4.1. The partial dataset resolved bbbb, resolved bbτ + τ − , bbγγ, and bbVV * analyses are described in Sections 4.2, 5.2, 6.1, and 7.1, respectively.
The limits for both benchmark models are shown in Figure 57.  Figure 57. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant ggF HH production combining channels using a partial Run 2 CMS dataset in the context of (a) spin-0 and (b) spin-2 benchmark models. The ±1σ and ±2σ uncertainties on the expected limits are shown with the colored bands. [79] Figure 57. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of ggF resonant HH production combining channels using a partial Run 2 CMS dataset in the context of (a,b). The ±1σ and ±2σ uncertainties on the expected limits are shown with the colored bands [79]. Figure 58 shows an overlay of the limits from the combination, the individual channels used in the combination, and the partial Run 2 bb νqq analysis [82].  Figure 58. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant ggF HH production using a partial Run 2 CMS dataset in the context of (a) spin-0 and (b) spin-2 benchmark models for a combination as well as individual channels. The notation (2+2) refers to resolved bb¯bb¯ reconstruction using two pairs of jets, (2+1) refers to semi-boosted reconstruction using one large-R jet and one pair of small-R jets and (1+1) refers to boosted reconstruction using two large-R jets. [ . Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant ggF HH production using a partial Run 2 CMS dataset in the context of (a,b) for a combination as well as individual channels. The notation (2+2) refers to resolved bbbb reconstruction using two pairs of jets, (2+1) refers to semiboosted reconstruction using one large-R jet and one pair of small-R jets, and (1+1) refers to boosted reconstruction using two large-R jets [83].

Searches for Resonant SH and SS Production
Both ATLAS and CMS have conducted searches for resonant production of additional scalar bosons. The searches assume ggF production of a heavy scalar that decays into a Higgs boson plus an additional scalar S (SH) or into a pair of S (SS). As discussed in Section 2.2, S is expected to have mass-dependent couplings to SM particles that are proportional to those of the SM Higgs boson, with a possible suppression of couplings to vector bosons. The branching ratios of SH and SS as a function of m S are shown in Figures 59 and 60, respectively, assuming no vector boson coupling suppression. Figure 58. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant ggF HH production using a partial Run 2 CMS dataset in the context of (??) spin-0 and (??) spin-2 benchmark models for a combination as well as individual channels. The notation (2+2) refers to resolved bbbb reconstruction using two pairs of jets, (2+1) refers to semi-boosted reconstruction using one large-R jet and one pair of small-R jets and (1+1) refers to boosted reconstruction using two large-R jets. [      The primary background process in the τ lep τ had final state is tt production with 1452 one real lepton and one real τ lepton, which is estimated using simulation. In the 1453 τ had τ had final state, the background is dominated by QCD multijet events with hadronic 1454 jets misidentified as τ had and Z+jets events where the Z decays into two real τ had .

1455
Background processes that contain two real τ had , either from the decay of a Z boson 1456 or two W bosons, are estimated using a τ-embedding method [85]. Data events that

CMS SH → bbτ + τ − Search
The CMS Collaboration conducted a search for resonant production of a Higgs boson plus an additional scalar in the bbτ + τ − channels in the bbτ had τ had and bbτ lep τ had final states [84] using the full Run 2 dataset of 137 fb −1 of data collected between 2016 and 2018. A heavy scalar H decays into a Higgs boson h and an another scalar h S with SM Higgs boson-like couplings to SM particles. The decay mode in which h decays to τ + τ − and h S decays to bb is used.
The search is interpreted in the context of the NMSSM as well as generic models predicting extra scalars. Limits are set for m H in the range 240 GeV to 3 TeV and m h S in the range 60 GeV to 2.8 TeV with m h S < m H − 125 GeV.
Events are selected with a combination of triggers that require a single lepton, a lepton plus a τ had-vis , or a pair of τ had-vis . Additionally, events are required to contain one lepton that passes isolation criteria and one τ had-vis that passes τ lepton identification criteria or two τ had-vis that pass the identification criteria. The τ lepton visible decay products are required to have opposite sign electric charge and be geometrically separated by ∆R > 0.5. Events are also required to contain at least one b-tagged anti-k t R = 0.4 jet plus an additional jet. If the event contains two or more b-tagged jets, the h candidate is built from the two that have the highest p T . If the event contains one b-tagged jet, it is paired with the remaining jet with the highest b-tagging score to reconstruct h.
A set of multiclassifier neural networks (NNs) are used to discriminate signal hypotheses from background. Each NN is trained on a group of signal hypotheses with similar m H and m h S values. The NNs use kinematic properties of the selected particles, invariant masses of combinations of particles, and the b-tagging scores of the jets used to reconstruct h. Each NN classifies events into five exclusive categories: signal events, background events with two real τ leptons, events with hadronic jets misidentified as τ had , tt events with no more than one real τ had , and events from other background processes such as Z+jets and single Higgs boson production.
The primary background process in the τ lep τ had final state is tt production with one real lepton and one real τ lepton, which is estimated using simulation. In the τ had τ had final state, the background is dominated by QCD multijet events with hadronic jets misidentified as τ had and Z+jets events where the Z decays into two real τ had . Background processes that contain two real τ had , either from the decay of a Z boson or two W bosons, are estimated using a τ-embedding method [85]. Data events that contain two muons are selected, and the muon signatures are replaced with signatures of simulated τ had with the same kinematic properties.
The QCD multijet background as well as tt and W+jets events that contain hadronic jets misidentified as τ had are estimated using a data-driven fake factor method. Three dedicated control regions (CRs) enriched in QCD multijet, W+jets and tt events are used to determine the fake factors associated with each process. The QCD multijet CR is defined by a requirement that the two τ lepton visible decay products have the same sign electric charge. A large transverse mass of the lepton and E miss T requirement and a veto on b-tagged jets are used to define the W+jets CR. Events in the tt CR are required to have more than two leptons. Corrections to the fake factors are derived in additional control regions to account for nonclosure effects.
All other background processes for both final states are estimated using simulation. Experimental and theoretical uncertainties are evaluated on the background estimates as well as the signal hypotheses. Uncertainties are evaluated on several aspects of the detector response modeling, the dominant of which is the τ had trigger efficiency and the rate at which leptons are misidentified as τ had . The τ-embedding process has associated uncertainties, including an uncertainty on the composition of the events. The fake factor method has uncertainties that include normalization and nonclosure effects.
An extended binned likelihood is used to fit the outputs of the NNs in each of the five categories. No significant excesses are observed over the predicted background. The upper limits on the production cross-section times branching ratio are set on the signal hypotheses at 95% CL using the CL s method, combining both final states. The limits are shown in Figure 61 as a function of m h S for each value of m H . The limits are also shown in the m H − m h S plane in Figure 62 along with an indication of the regions that exclude the maximally allowed values in the NMSSM.   Figure 61. Expected and observed 95% confidence level limits on the production cross-section times branching ratio of resonant production of a Higgs boson plus an additional scalar in the bbτ + τ − channel using the full Run 2 CMS dataset. The limits are shown as a function of the mass of the additional scalar m h S and are scaled by orders of ten to show the different masses of the heavy scalar m H . The ±1σ and ±2σ uncertainties on the combined expected limits are shown with the colored bands [84].

ATLAS SS → WW ( * ) WW ( * ) Search
For extended Higgs sector models in which bosonic decay modes of S are not suppressed, S → WW ( * ) is the dominant decay mode for m S 135 GeV. Therefore, the WW ( * ) WW ( * ) channel provides a combination of a clean signature with low background and a high branching ratio to search for resonant X → SS production.
The search for HH in the WW * WW * channel conducted by the ATLAS Collaboration using 36.1 fb −1 data as described in Section 8.1 was also used to search for SS production in the same channel. The methodology is the same for both signal models, with dedicated selection optimizations for each X → SS signal hypothesis. The event yields in each category after the selection optimized for the m X = 280 GeV m S = 135 GeV signal hypothesis are shown in Figure 63. The upper limits on the production cross-section times branching ratio are set at 95% CL using the CL s method on signal hypotheses with m X = 340 GeV and m S in the range of 135-165 GeV as well as m S = 135 GeV and m X in the range of 280-340 GeV. The limits are shown in Figure 64.   The CMS analysis uses the bbτ + τ − final state to search for SH production where H decays to τ + τ − and S decays to bb. The search is performed over a large range of m S from 60 to 2975 GeV. For m S < m H , the branching ratio of S decaying to bb remains approximately the same as for HH. Under models that follow the sum rule, S → bb remains important for m S > m H . Constraining the search to use S → bb, reduces the total branching ratio by a factor of about 2 but offers advantages over including S → τ + τ − decays. For low values of m X and m S , an S → τ + τ − decay would be very soft, making it difficult to select events with triggers that require leptons or hadronically decaying τ leptons. Additionally, the m H constraint allows the τ + τ − system to be approximately reconstructed as a way to discriminate signal from background. The exclusion limits are presented in a two-dimensional mass plane under the assumption of the NMSSM.
The ATLAS analysis uses the WW * WW * decay mode in multilepton final states to search for SS production. The search is constrained to m S between 135 and 165 GeV. Under the assumption that the sum rule is not followed, S → WW * has the dominant branching ratio for m S > 135 GeV. The upper bound is used based on the assumption that for m S > 2m t , S would decay primarily into tt. This assumption has been shown to be invalid. Future iterations of the search would benefit from an extended mass range. Additionally, the exclusion limits in future iterations would be more useful if presented in a two-dimensional plane.
These are the first two experimental results testing the existence of such models. A wide range of final states and topologies need to be studied to provide a comprehensive search for both SH and SS production. The two searches represent the beginning of a much larger search program for both collaborations.

Conclusions
A comprehensive overview of the state-of-the-art experimental searches for resonant production of scalar boson pairs by the ATLAS and CMS Collaborations using √ s = 13 TeV pp collisions from Run 2 LHC data is presented. The searches use between 35.9 and 139 fb −1 of data. Searches for resonant HH production are performed in numerous channels by both ATLAS and CMS. The results from the searches are interpreted in the context of multiple benchmark signal models. A summary of the searches performed by ATLAS is given in Table 3 and a summary of the searches performed by CMS is given in Table 4. Results from the individual channels are combined and shown in Sections 9.1 and 9.2.  Additionally, a search for SH production in the bbτ + τ − channel is performed by CMS, and a search for SS production in the WW * WW * channel is performed by ATLAS.
No significant excesses have been observed, but searching for resonant scalar pair production is a very active area of research with many ongoing analyses in both ATLAS and CMS. Many updated and new search results are expected as the full Run 2 datasets continues to be analyzed and as the LHC Run 3 data collection begins.