Combination of searches for resonant Higgs boson pair production using $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector

A combination of searches for a new resonance decaying into a Higgs boson pair is presented, using up to 139 fb$^{-1}$ of $pp$ collision data at $\sqrt{s}=13$ TeV recorded with the ATLAS detector at the LHC. The combination includes searches performed in three decay channels: $b\bar{b}b\bar{b}$, $b\bar{b}\tau^+\tau^-$ and $b\bar{b}\gamma\gamma$. No excess above the expected Standard Model background is observed and upper limits are set at the 95% confidence level on the production cross section of Higgs boson pairs originating from the decay of a narrow scalar resonance with mass in the range 251 GeV-5 TeV. The observed (expected) limits are in the range 0.96-600 fb (1.2-390 fb). The limits are interpreted in the Type-I Two-Higgs-Doublet Model and the Minimimal Supersymmetric Standard Model, and constrain parameter space not previously excluded by other searches.

The discovery of a Higgs boson (ℎ) with a mass of about 125 GeV at the Large Hadron Collider (LHC) in 2012 [1,2] has initiated a major effort to understand the possible connections of this particle with new physics phenomena.In particular, the possibility that heavy particles beyond the Standard Model (SM) may couple to the Higgs boson has attracted much attention.This possibility is particularly compelling because scalar particles like the Higgs boson are dimension-one operators.Such operators allow couplings between new particles and the Higgs boson that are unsuppressed by higher energy scales [3].These couplings could give rise to sizable effects evading indirect constraints from other searches.Related signatures in which a heavy particle decays into a Higgs boson pair (ℎℎ) are commonplace in almost all extensions of the Higgs sector, including models with additional weak isospin singlets, doublets or triplets [4][5][6][7], and in more exotic models, e.g.models with extra dimensions [8].This type of resonant ℎℎ production typically has a cross section that is several orders of magnitude larger than the SM non-resonant ℎℎ production [9], depending on the model.This Letter presents the combination of searches for resonant Higgs boson pair production in three ℎℎ decay channels: four -jets ( b b), two -jets and two -leptons ( b +  − ), and two -jets and two photons ( b).These decay modes provide better sensitivity to Higgs boson pair production than other channels because they either have large Higgs boson branching ratios or are easy to discriminate from background processes.The searches under consideration use the LHC Run-2 dataset corresponding to an integrated luminosity of up to 139 fb −1 of proton-proton ( ) collision data recorded with the ATLAS detector at a centre-of-mass energy of √  = 13 TeV.The ATLAS experiment at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and nearly 4 coverage in solid angle [10][11][12], and uses an extensive software suite [13] for both simulation and data analysis.
Previous combinations of searches for resonant Higgs boson pair production were performed using only part of the LHC Run-2 dataset corresponding to up to 36.1 fb −1 of   collision data.The ATLAS Collaboration combined the decay modes  b b,  b +  − ,  b +  − ,  +  −  +  − ,  b and  +  −  [14], and the CMS Collaboration combined the decay modes  b b,  b +  − ,  b ,  b +  − and  b [15].More recently, new searches were performed using the LHC Run-2 dataset in the  b b,  b +  − and  b decay channels by the ATLAS Collaboration [16][17][18] and in the  b b,  b +  − ,  b +  − ,  +  −  +  − ,  +  −  +  − and  +  −  +  − decay channels by the CMS Collaboration [19][20][21].These searches focus on a particle with mass greater than 250 GeV produced via gluon-gluon fusion whose natural width is small compared with the experimental mass resolution.Gluon-gluon fusion is chosen because it is typically a dominant production mechanism for this signature in most extensions of the Higgs sector.
Monte Carlo simulated signal samples for a new, narrow scalar resonance  that is produced via gluon-gluon fusion and decays into ℎℎ ( → ℎℎ), are generated at leading order accuracy in the strong coupling constant using the MADGRAPH5_aMC@NLO v 2.6.1 package [22] with the NNPDF2.3loparton distribution function set [23].The event generator is interfaced with the Herwig v7.1.3program [24,25] to model the parton shower, hadronisation and underlying event.Different samples are generated in which the  particle mass,   , is varied in the range 251 GeV-5 TeV.Its natural width is set to 10 MeV, which is negligible compared with the experimental mass resolution, and the mass of the Higgs boson,  ℎ , is set to 125 GeV.The interference with non-resonant ℎℎ production is neglected.The interference is negligible for the combined cross-section upper limits that are presented in this Letter, due to the narrowness of the assumed  resonance width.For the interpretation, neglecting the interference is an approximation.Descriptions of the data and simulation samples used to model backgrounds from SM processes can be found in Refs.[16][17][18].
Each of the three decay channels included in this combination provides the highest sensitivity in the search for resonant ℎℎ production for a portion of the   range considered.The  b b final state has the highest branching ratio (33.9%) and is the most sensitive channel for high resonance masses (  ≳ 800 GeV), but the need to impose high  T -jet trigger thresholds, the ambiguity in the assignment of -jets to Higgs boson candidates, and the large multĳet background reduce its sensitivity at lower   .The  b +  − decay has a moderate branching ratio (7.3%) and moderate background contamination and provides the best sensitivity in the intermediate resonance mass range (350 GeV ≲   ≲ 800 GeV).The  b final state has the lowest branching ratio (0.3%), but benefits from the high trigger and reconstruction efficiencies for photons.Combined with an excellent diphoton mass resolution, this channel has the best sensitivity at low resonance masses (  ≲ 350 GeV).All three channels require the presence of -quark initiated jets (-jets), which are identified using the DL1r algorithm [26,27] at the 77% efficiency working point.The DL1r algorithm is a deep-neural-network classifier that uses track and vertex information to classify the jet with respect to whether it contains -or -hadrons.At the 77% -jet identification efficiency operating point, light-jet (-jet) rejection factors of 170 (5) are achieved in a sample of simulated SM  t events.
The  b b analysis [16] is performed in two complementary channels: the resolved channel in which Higgs boson candidates are formed from four small-radius ( = 0.4) anti-  jets [28,29], and the boosted channel in which the high- T Higgs boson candidates ( T > 250 GeV) are reconstructed as two separate large-radius ( = 1.0) anti-  jets.The resolved channel relies on a combination of small-radius -jet triggers [30,31] and uses 126 fb −1 of data, while the boosted channel uses large-radius jet triggers [30,32] without any -tagging requirement applied and uses 139 fb −1 of data.In the resolved channel, events are required to have at least four  = 0.4 anti-  -jets [33,34], and a boosted decision tree (BDT) is used to pair the -jets to form the two Higgs boson candidates.In the boosted channel events are placed in the 2, 3 and 4 categories according to the number of -tagged track-jets associated with each of the two  = 1.0 anti-  jets that form the Higgs boson candidates.These track-jets are built from inner-detector tracks clustered using the anti-  algorithm with a variable radius and are matched to the large-radius jets using the ghost association method [35], as described in Ref. [16].The invariant masses of the two Higgs boson candidates are used to classify events into signal, validation, and control regions.The final observable in both channels is provided by the invariant mass of the two-Higgs-boson system, which is reconstructed with a resolution of ∼ 5-6%, thanks to the application in the analysis of a kinematic rescaling of the ℎ →  b system to the ℎ-boson mass that improves the reconstructed ℎℎ mass resolution.The search is performed in the range 251 GeV <   < 1.5 TeV (900 GeV <   < 5 TeV) in the resolved (boosted) channel.The two channels are made orthogonal by explicitly vetoing resolved events that appear in the boosted selection, and are statistically combined in the overlapping mass region to enhance the sensitivity of the search.
The  b +  − analysis [17] defines different event categories depending on decay modes of -leptons.The -leptons decay into final states containing hadrons ( had ) or leptons ( lep ) [36,37], and the analysis considers Higgs boson decays into  had  had and  lep  had .The  had  had channel relies on single- had and di- had triggers [30,32], while the  lep  had channel uses single-lepton and lepton-plus- had triggers [30,32,38,39].The search exploits the LHC Run-2 dataset of 139 fb −1 .Events in the  lep  had channel are categorised in two orthogonal regions depending on whether they satisfy a single-lepton trigger or a lepton-plus- had trigger.In the  had  had category events are required to have exactly two reconstructed  had candidates with opposite charge and no reconstructed electrons or muons, while in the  lep  had category events are required to have exactly one electron or muon and one  had with opposite charge.In both categories exactly two  = 0.4 anti-  -jets are required.The analysis uses the output score of a mass-parameterised neural network [40] as the final observable in each signal region category.The neural network can discriminate the mass of a possible signal excess with a resolution of 5-10%, depending on   .All  had  had and  lep  had signal regions are statistically combined in the search for resonances in the range 251 GeV <   < 1.6 TeV.
The  b analysis [18] relies on single-photon and diphoton triggers [30,39] and uses the LHC Run-2 dataset of 139 fb −1 .Events are required to have exactly two photons [41], two  = 0.4 anti-  -jets, and zero electrons or muons.Two BDTs are trained to separate the signal from backgrounds: one is trained against the  t background and another is trained against the single Higgs boson background.The signal region selection requires the mass of the diphoton system to be compatible with the SM Higgs boson mass and the mass of the  b system to be compatible with the expectation for the corresponding signal events.In addition, different requirements are applied on a combination of the two BDT scores with a tuned coefficient across the mass hypotheses.For this combination, the  b search follows the analysis strategy detailed in Ref. [18] but extends the 251 GeV <   < 1.0 TeV mass range used there to 251 GeV <   < 1.3 TeV.The final discriminating variable in this channel is the diphoton invariant mass, which is reconstructed with a resolution of about 1%.The resolution of the mass of the ℎℎ system in this channel is ∼ 2-3%, for the mass range considered in this search, thanks to the excellent di-photon mass reconstruction and the use of a rescaling of the reconstructed ℎℎ mass.
The results of the combination presented in this Letter are obtained from a likelihood function (, ì ), where  represents the parameter of interest of the model, namely the resonant ℎℎ production cross section, and ì is a set of nuisance parameters, including the systematic uncertainty contributions and background parameters that are constrained by control regions in data.The global likelihood function  (, ì ) is obtained as the product of the likelihoods of the three searches included in the combination, which are themselves products of likelihoods computed from the final observables in the single analysis categories.In this combination, the  b b,  b +  − and  b analysis signal regions are either orthogonal due to the different number and type of physics objects required in their final states, or have a negligible number of overlapping events.They are therefore treated as statistically independent.The profile-likelihood-ratio test statistic [42] is used to obtain upper limits on the cross section of resonant ℎℎ production with the CL s method [43].For signal masses up to 3 TeV, the limits are computed using asymptotic formulae [42].At higher masses, the asymptotic approximation is inaccurate due to the limited number of background candidates and the limits are computed by sampling pseudo-experiments.
Systematic uncertainties related to the data-taking conditions, such as those associated with the integrated luminosity and the modelling of the effect of multiple interactions in the same and neighbouring bunch crossings (pileup), are considered correlated across the searches.Uncertainties related to reconstructed quantities used by multiple searches and theoretical uncertainties on simulated signal and background processes are treated as correlated where appropriate.Part of the -tagging uncertainties and the signal parton shower modelling uncertainties are treated as uncorrelated due to different definitions and implementations in the individual analyses.Different uncertainties dominate the combined limit in different   ranges.These correspond to the dominant uncertainties of the channel most sensitive to resonant masses in that range.They are: photon scale and resolution, signal and background modelling from  b for low masses [18];  plus heavy-flavour jets modelling and single-Higgs boson plus heavy-flavour jets modelling from  b +  − for intermediate masses [17]; and multĳet modelling, jet mass resolution and flavour-tagging from  b b for high masses [16].The total impact of the systematic uncertainties on the combined limit is between 3% and 10%, depending on   , and the impact of ignoring systematic uncertainty correlations between the three channels is of the order of 1%.
The results of the combination are compatible with the predictions from the SM backgrounds across the tested   range from 251 GeV to 5 TeV. Figure 1(a) shows the local -value as a function of   for a narrow resonance that decays into a pair of SM Higgs bosons.The largest deviation is observed at 1.1 TeV and corresponds to a local significance of 3.3 , which is driven mainly by the  b +  − channel.The global significance of this excess is estimated to be 2.1  using the method outlined in Ref. [44]. Figure 1(b) shows the upper limits at the 95% confidence level (CL) on the resonant ℎℎ production cross section as a function of   , assuming that ℎ is the SM Higgs boson.The observed (expected) upper limits are in the range 0.96-600 fb (1.2-390 fb), depending on   .The  b search is the most sensitive at low mass, the  b +  − search is the most sensitive in the 350-800 GeV range, and the  b b search dominates for high resonant masses, demonstrating the complementarity of these three searches.These results represent an improvement in the upper limits by a factor between two and five, depending on   , relative to the previous ATLAS combination [14].The improvement comes not only due to the use of a larger dataset but also due to improvements in object identification and analysis techniques in background estimation.The upper limits are interpreted in the context of the Type-I Two-Higgs-Doublet Model (2HDM) [5,45] and the Minimal Supersymmetric SM (MSSM) [46][47][48][49][50].For these interpretations, the branching ratios of the Higgs boson are fixed to the predictions of the model under consideration.
The 2HDM is an extension of the SM with an additional Higgs doublet, leading to five Higgs bosons after electroweak symmetry breaking, three neutral and two charged.The neutral Higgs bosons, assuming charge-parity (CP) conservation, can be CP-even (ℎ and , with   >  ℎ ) or CP-odd ().The resonant ℎℎ production is assumed to come from  → ℎℎ for this interpretation.The ℎ boson is assumed to be the Higgs boson observed at the LHC with  ℎ = 125 GeV.The other Higgs bosons (, ,  ± ) are assumed to be mass degenerate, i.e.,   =   =  ± , and the Higgs potential parameter  2  12 is fixed to , where tan  is the ratio of the vacuum expectation values of the two Higgs doublets.No constraint is imposed on cos( − ), where  is the mixing angle between the CP-even Higgs bosons.It, along with tan  and   , is taken as a free parameter.Widths and branching ratios are calculated using the 2HDMC program [51].The  boson gluon-gluon fusion production cross section calculation uses the SusHi package [52,53], which includes corrections up to next-to-next-to-leading order in   [54][55][56], massive quarks [57,58] and electroweak corrections by light fermions [59,60].The procedure for the theoretical calculations follows Ref. [9].The upper limits on the Higgs boson pair production cross section are interpreted as constraints on two benchmark planes:   -tan  for given cos( − ) values, shown in Figures 2(a)-2(b), and cos( − )-  for fixed tan  values, shown in Figures 2(c)-2(d).The interpretation is given in the Type-I 2HDM in which gluon-gluon fusion is the dominant production mechanism throughout the parameter space where the searches have sensitivity.The  boson has a finite natural width in the 2HDM.This is taken into account by considering in the combination only channels for which the experimental mass resolution is large in comparison to the  boson natural width.In particular, the  b ( b b and  b +  − ) channel upper limits are valid for  boson natural widths Γ  /  up to 2% (5%).The results are sensitive to cos( − ) values that are not probed by the SM Higgs boson coupling measurements [61].For example, the point tan  = 10, cos( − ) = −0.1 is excluded at 95% CL for   values in the range 270-810 GeV, a region allowed by Higgs boson coupling measurements.The  b channel limits are valid for  boson natural widths Γ  /  < 2%, so limits from this channel are quoted for the parameter space that satisfies this requirement.The limits in the  b b and  b +  − channels are only valid for Γ  /m  < 5%, so no limits are quoted in the dark shaded regions, where the width exceeds this value.
The MSSM has a Type-II 2HDM Higgs sector structure [62] and, therefore, includes the 2HDM parameters discussed earlier.The resonant ℎℎ production is assumed to come from  → ℎℎ for this interpretation.
Supersymmetry constrains the number of free parameters in the Higgs sector at lowest order to be two, taken here as   and tan .Radiative corrections have a large impact on the MSSM and the limits are influenced by how the supersymmetry parameters are chosen, with each choice defining a particular scenario.In this Letter, the M 125 ℎ,EFT and M 125 ℎ,EFT ( χ) scenarios are used [63].These scenarios have supersymmetry mass parameters that are not related to the Higgs sector at a very high energy scale such that the low tan  region has  ℎ ≈ 125 GeV.The M 125 ℎ,EFT ( χ) scenario differs from M 125 ℎ,EFT in that it includes low-mass neutralinos and charginos.The MSSM gluon-gluon fusion cross section is obtained using a SusHi implementation that includes supersymmetric QCD corrections [64,65].The Higgs boson masses and mixing are calculated with the Feynhiggs program [66][67][68][69][70][71][72][73] as described in Ref. [74].The Higgs boson branching ratio calculation uses FeynHiggs, HDECAY and PROPHECY4f as described in Refs.[75,76].The upper limits on the Higgs boson pair production cross section are interpreted in the MSSM as constraints in the   -tan  plane, as shown in Figure 3.The combined analysis excludes parameter space in the region 2 ≲ tan  ≲ 5, which is not excluded by the / → ,  →  ℎ,  →   or  ± →  searches [77][78][79][80].
Figure 3: Observed and expected exclusion limits at the 95% CL on the MSSM parameter space for the (a) M 125 ℎ,EFT and (b) M 125 ℎ,EFT ( χ) benchmark scenarios for each of the individual channels and their combination.The  b channel limits are valid for  boson natural widths Γ  /  < 2%, so limits from this channel are quoted for the parameter space that satisfies this requirement.The limits do not apply in the dark shaded regions, where the mass of the lightest CP-even Higgs boson, ℎ, is not compatible with 125 GeV within the experimental mass resolution.
In conclusion, this Letter reports a combined search for a narrow-width resonance decaying to ℎℎ in the  b b,  b +  − and  b final states using up to 139 fb −1 of   collision data at √  = 13 TeV collected by the ATLAS experiment.The data are found to be consistent with the SM prediction and upper limits at the 95% CL are set on the resonant  → ℎℎ cross section assuming   to be in the range 251 GeV-5 TeV.The observed (expected) upper limits are in the range 0.96-600 fb (1.2-390 fb) and they constitute an improvement of a factor of 2-5, depending on   , with respect to the previous ATLAS combined result [14].The results are also interpreted in the context of the Type-I 2HDM and MSSM, excluding parameter space that was hitherto allowed by the most sensitive searches for these models.

Figure 1 :
Figure 1: (a) Local -value and (b) observed and expected upper limits at the 95% CL on the resonant Higgs boson pair production cross section as a function of the resonance mass   .The symbol ℎ denotes a SM Higgs boson with a mass 125 GeV.

Figure 2 :
Figure 2: Exclusion limits at the 95% CL on the Type-I 2HDM parameter space.Observed and expected limits for the combination of all channels and expected limits for each of the individual channels are presented.The 2HDM parameters are shown in the   -tan  plane for (a) cos( − ) = −0.1 and (b) cos( − ) = 0.1, and in the cos( − )-  plane for (c) tan  = 1 and (d) tan  = 10.The  b channel limits are valid for  boson natural widths Γ  /  < 2%, so limits from this channel are quoted for the parameter space that satisfies this requirement.The limits in the  b b and  b +  − channels are only valid for Γ  /m  < 5%, so no limits are quoted in the dark shaded regions, where the width exceeds this value.