Search for resonant pair production of Higgs bosons decaying to bottom quark-antiquark pairs in proton-proton collisions at 13 TeV

A search for a narrow-width resonance decaying into two Higgs bosons, each decaying into a bottom quark-antiquark pair, is presented. The search is performed using proton-proton collision data corresponding to an integrated luminosity of 35.9 fb−1 at s=13\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s}=13 $$\end{document} TeV recorded by the CMS detector at the LHC. No evidence for such a signal is observed. Upper limits are set on the product of the production cross section for the resonance and the branching fraction for the selected decay mode in the resonance mass range from 260 to 1200 GeV.


Introduction
The discovery of a Higgs boson (H) [1-3], with mass of 125 GeV [4,5] and properties consistent with the standard model (SM) of particle physics at the CERN LHC, motivates searches for resonances via their decays into Higgs bosons. Several theories for physics beyond the SM posit narrow-width resonances decaying into pairs of Higgs bosons (HH). For instance, models with a warped extra dimension [6] predict the existence of new particles such as the spin-0 radion [7-9] and the spin-2 first Kaluza-Klein (KK) excitation of the graviton [10-12], which could decay to HH. These models have an extra warped spatial dimension compactified between two branes, with an exponential metric κl, κ being the curvature and l the coordinate of the extra spatial dimension [13]. The benchmark of the model is the ultraviolet cutoff of the theory, Λ ≡ √ 8πe −κl M Pl , M Pl being the Planck scale. In proton-proton (pp) collisions at the LHC, the graviton and the radion are produced primarily through gluon-gluon fusion and are predicted to decay to HH with branching fractions of approximately 10 and 23%, respectively [14].
Previous searches for resonant HH production have been performed by the ATLAS and CMS Collaborations with pp collisions at √ s = 8 and 13 TeV. The decay channels studied include bbbb [15][16][17], bbττ [18], bbγγ [19,20], γγWW [21], and bbWW [22]. This paper reports the results of a search for narrow-width resonances in the 260-1200 GeV mass range, decaying into a pair of Higgs bosons, each decaying into a pair of bottom quarks. The search is performed using pp collision data collected at √ s = 13 TeV with the CMS detector at the CERN LHC, corresponding to an integrated luminosity of 35.9 fb −1 . The main challenge of this search is to discriminate the final-state signature of four bottom quark jets from the overwhelming multijet quantum chromodynamics (QCD) background. This is addressed by dedicated online selection criteria that include b jet identification and by a model of the multijet background that is tested in control regions of data. The analysis closely follows the approach adopted for the 8 TeV data [15] but the sensitivity for high resonance mass values is enhanced because of the significant increase in production cross section at 13 TeV, a new trigger strategy and a more efficient algorithm for identifying jets originating from bottom quarks.

Detector and simulated samples
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. A silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections, reside within the solenoid. Forward calorimeters extend the pseudorapidity (η) [23] coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [23] Simulated samples of signal events are produced using various Monte Carlo (MC) event generators, with the CMS detector response modeled using the GEANT4 [24] program. To model the production of a generic narrow-width spin-0 resonance, we use an MC simulation of the bulk radion produced through gluon fusion. The momenta and angular distributions for decay products of a spin-2 resonance are distinct from those for a spin-0 resonance, and result in different kinematic distributions. Therefore, we evaluate the signal efficiencies for a narrowwidth spin-2 resonance from a separate simulation of the first excitation of a bulk KK graviton produced through gluon fusion and forced to decay to a pair of Higgs bosons with the parame-optimizing for the expected sensitivity and takes into account the uncertainties associated with the background modeling. The mass range above 1200 GeV (high-mass region) is not covered by this search. Above 900 GeV, the Higgs bosons have a momentum considerably higher than their mass and the Higgs to bb decays are reconstructed more efficiently as single hadronic jets with a larger anti-k T distance parameter (0.8) [38].
Events are selected online by combining two different trigger selections to identify b jets, both using the CSV algorithm. For the first trigger selection, four jets with p T > 30 GeV and |η| < 2.4 are required. The latter requirement ensures that the jet lies within the tracker acceptance. Of those four jets, two are required to have p T > 90 GeV and at least three jets are required to be tagged as b jets. The second trigger selection requires four jets with p T > 45 GeV and at least three of those jets identified as b jets.
Events are selected offline by requiring at least four b tagged jets with p T > 30 GeV and |η| < 2.4. The selected jets are combined randomly into pairs to form two Higgs boson candidates with masses m H 1 and m H 2 . For the LMR, HH candidates are chosen from the four selected jets such that |m H − 120 GeV| < 40 GeV for each candidate Higgs boson. For the MMR the H candidates are selected using ∆R = √ (∆η) 2 + (∆φ) 2 less than 1.5, where ∆η and ∆φ are the differences in the pseudorapidities and azimuthal angles (in radians) of the two jets.
In the two-dimensional space defined by the reconstructed masses of the two Higgs boson candidates, H 1 and H 2 , a circular signal region (SR) is defined with R < 1, where R is defined as: The central mass value (M) is the average of the means of the m H 1 and m H 2 distributions for simulated signal events and the parameter r is set to 20 GeV. The centers of these circular regions have been determined separately for the LMR and MMR and found to be 120 and 125 GeV, respectively. If there are multiple HH candidates in an event, the combination that minimizes R 2 is used.
After these event selection criteria are applied, the dijet invariant mass resolution for m H 1 and m H 2 is approximately 10-13%, depending on the p T of the reconstructed Higgs boson, with a few percent shift in the value of the mass peak, relative to 125 GeV. The Higgs boson mass resolution is further improved by applying multivariate regression techniques similar to those used in the searches for SM Higgs bosons decaying to bb in CMS [39,40]. The regression estimates a correction that is applied after the standard CMS jet energy corrections [34,41], and it is computed for individual b jets to improve the accuracy of the measured energy with respect to the b quark energy. To this end, a specialized boosted decision tree [42] is trained on simulated b jets from tt events, with inputs that include observables related to the jet structure and b tagging information. The average improvement in the Higgs boson mass resolution, measured with simulated signal samples, is 6-12%, depending on the p T of the reconstructed Higgs boson. The use of the regression technique increases the sensitivity of the analysis by 5-20% depending on the mass hypothesis. The regression technique is validated with data samples of Z → (ee, µµ) events with two b tagged jets and in tt-enriched samples [39]. The cumulative selection efficiencies of the selection criteria described above for the graviton and radion signal benchmarks are reported in Fig. 1.
The reconstructed resonance mass (m X ), computed as the invariant mass of H 1 and H 2 , is displayed for simulated signal events with different mass hypotheses in Fig. 2. In order to improve the resolution for the resonance mass, the momenta of the reconstructed b quarks are corrected by constraining the invariant mass of the Higgs boson candidates to be 125 GeV. Since jet direction is reconstructed with better resolution than jet p T , this constraint mainly benefits the latter. The improvement in resolution for the reconstructed signal resonance ranges from 20 to 40% depending on the mass hypothesis, resulting in an improvement of the sensitivity by 10-20%. The shift in the reconstructed m X , due to these corrections is also shown in Fig. 2. The mass shift is linear in m X and is caused by the asymmetry of the corrections due to the jet momentum resolution across the p T range considered.

Signal and background modeling
To search for signal events with various mass hypotheses we fit the m X distribution for data events in the SR to the sum of two parametric models, one for the signal and one for the SM background, which is dominated by multijet production. This procedure is performed independently in three regions: two within the LMR, as described below, and the MMR. Quantum interference between the signal and non-resonant SM background processes is neglected.
The parametric signal models are built by fitting the m X distributions obtained from simulated signal events. The shape of the signal m X distribution is different for the LMR and the MMR, as a result of the different event selection criteria. In the LMR, a sum of two Gaussian functions is used, to account for the effect of reconstructing the two Higgs boson candidates using jets not originating from their decays. This occurs in about 30% of the events for the lowest resonance mass hypothesis, decreasing to about 5% at the highest resonance mass hypothesis. Five parameters are required, the mean and width of the two Gaussian functions and the ratio of their integrals. In the MMR, the signal is modeled with a function that has a Gaussian core smoothly extended on both sides via exponential tails. This requires two parameters for the mean and width of the Gaussian function, and two parameters for the exponential tails [15].
The parametric function used to model the background m X distribution is obtained from control regions in data in which events are expected to have kinematic properties similar to events in the SR. These control regions are: (i) sideband regions (SB) defined as 1 < R < 2 and  Fig. 3), and (ii) events in the SR and both SB regions that do not qualify as signal events due to a requirement that one of the four jets used to reconstruct m X is not identified as a b jet (anti-tag selection). For the LMR selection, the m X distributions in all control regions exhibit a kinematic turn-on, due to the trigger requirements, followed by a smoothly falling distribution at the larger m X values. The m X distributions in all control regions for the MMR selection exhibit a similar smoothly falling shape with increasing m X . This feature allows the adoption of a common form of the parametric background model in the SR and control regions. The parameters of the background model are determined by the fit used to extract a possible signal. In the following we describe the derivation of the background model from the different control regions. For the LMR, it is difficult to model simultaneously the turn-on and the tail of the m X distribution. Thus, the modeling of the background is split into two ranges, LMR1 and LMR2, with two different parametric models to accommodate either a turn-on shape or a falling distribution, respectively.
An indirect method of determining the boundary between LMR1 and LMR2 is used to avoid revealing the m X distribution in the SR before fitting for a possible signal. The method uses two selections in the SR and SB, as described in Table 1. For the background, the shape of the m X distribution in the SR (A) is predicted from the shape of the m X distribution in the SR using the anti-tag selection (C) but reshaped by the ratio, as a function of m X , between the m X distribution in the SB with the nominal selection (B) and the m X distribution in the SB with the anti-tag selection (D), so that A = CB/D. To test the validity of this method, new SR and SB regions, centered at 150 GeV, are defined along the diagonal in Fig. 3 in a similar fashion to those centered at 120 GeV. In that case the prediction E = GF/H for the m X distribution in the SR region centered at 150 GeV can be tested by directly comparing to the m X distribution in that region. Figure 4 shows the agreement between this prediction and the observed distribution, validating the use of this method for the SR distribution centered at 120 GeV. The boundary between LMR1 and LMR2 is then set at m X = 310 GeV, by optimizing for the expected sensitivity. , while the function defined in Ref. [43] Eq. (9) is used in LMR2 and the MMR. This function was originally used to describe a Compton spectrum and has three free parameters describing the mean, width and extent of the right side tail. In each case, the goodness of the fit, characterized by the χ 2 per degree of freedom, is found to be reasonable. The fit of the background model on the SB for the MMR is shown in Fig 5. In the absence of a theoretical prediction for the background model, an alternative one based on the Crystal Ball function [44] has been studied, and the difference between them drives the assessment of a systematic uncertainty due to the choice of the particular model. Pseudodatasets are generated from the alternative function and fitted with the nominal functions to compute biases in the reconstructed signal strength. This procedure is performed for each mass hypothesis and the corresponding biases range between 30-80 fb for the LMR, and 0.1-4 fb for the MMR.
The SM top quark pair production in the SR is estimated from simulation to contribute up to 10 and 15% of the selected events in the LMR and the MMR, respectively. Since the tt contribution to the background exhibits a shape very similar to that for the multijet process, it is implicitly included in the data driven estimate.  Figure 4: The predicted m X distribution in data for the LMR (squares) and the observed distribution in a SR (circles) centered at a Higgs boson mass of 150 GeV.

Systematic uncertainties
Sources of systematic uncertainties that affect the signal yields are listed in Table 2. The signal yield for a given production cross section is affected by a 2.5% systematic uncertainty in the measurement of the integrated luminosity at CMS [45]. The uncertainty in the signal normalization caused by the choice of the PDF set [46] contributes up to 3.5%. The jet energy scale [34,41] is varied within one standard deviation as a function of jet p T and η, and the efficiency of the selection criteria recomputed. The signal efficiencies are found to vary by up to 2.9%. The effect of the uncertainty in the jet energy resolution [34,41] is evaluated by smearing the jet energies according to the measured uncertainty. This causes variations in the signal efficiency of between 0.9 and 2.1%. The trigger efficiency is evaluated in a tt-enriched control sample in which an isolated muon is required in addition to the four b tagged jets with p T > 30 GeV and |η| < 2.4 demanded in the analysis. The efficiency of each online kinematic and b tagging requirement is evaluated separately and then the efficiency of all the trigger selection criteria is computed in terms of conditional probability. The associated systematic uncertainties are found to impact signal efficiencies at the level of 5-9%. Signal yields are corrected to match the b tagging efficiency measured in data [35]. The associated uncertainty is evaluated to be about 6-8%.
The systematic uncertainty associated with the choice of the parametric background model is evaluated as described in Section 5. We account for the bias as a signal-shaped systematic uncertainty in the background model with normalization centered at zero and a Gaussian uncertainty with standard deviation equal to the bias. The impact on the expected limit ranges between 0.3-1.5%.

Results
The m X distribution for data in the SR and results of the fit to the parametric background model are shown in Figs. 6, for the LMR1 and the LMR2, and 7 for the MMR. Parameters controlling

Source of
Impact in LMR (%) Impact in MMR (%) systematic uncertainty Signal Signal Luminosity 2.5 2.5 Jet energy scale 0.2-1.8 0.9-2.9 Jet energy resolution 0.9-2.1 1.0-1.5 b tagging scale factor 6.5-6.9 6.9-8.6 Trigger efficiency 6.4-9.0 5.3-7.0 PDF 1.5-2.2 2.1-3.5   Figure 9: The observed and expected upper limits on the cross section for a spin-0 resonance X → H(bb)H(bb) at 95% CL, using the asymptotic CL s method. The theoretical cross section for the production of a radion, with Λ = 3 TeV, κl = 35, and no radion-Higgs boson mixing, decaying to four b jets via Higgs bosons is overlaid. The transition between the LMR and the MMR is based on the expected sensitivity, resulting in the observed discontinuity.
the shapes and yields of the signal are allowed to float within ranges determined by the systematic uncertainties. The parameters and normalization of the multijet background shape are left to float freely in the fit. The shapes of the background-only fit are found to adequately describe the data in each of the search region.
The observed and expected upper limits on the cross section times branching fraction for pp → X → H(bb)H(bb) at 95% confidence level (CL) are computed using the modified frequentist CL s method [47][48][49][50]. These limits are shown in Figs. 8 and 9 for spin-2 and spin-0 hypotheses, respectively. The green and yellow bands represent the ±1 and ±2 standard deviation confidence intervals around the expected limits. The observed limits in these figures exhibit several deviations from the expected one beyond ±2 standard deviation. The local p-value of the most significant excess (deficit) is 2.6 (−3.6) standard deviation. In order to estimate the global significance of these deviations within the search range, we estimated, based on pseudo-experiments, a global probability to see three (four) positive or negative deviations from the expected limits consistent with narrow positive or negative signal contributions with the square sum of the significances exceeding 5.1 (5.6) standard deviation, as observed in data. In both cases, the global probability was found to be in a percent range, consistent with a statistical nature of the observed deviations. In Fig. 8, the NLO theoretical cross section for the gluon fusion production of a bulk KK-graviton with κl and κ set to 35 and 0.5M Pl respectively is shown. This KK-graviton is excluded at 95% CL in the mass ranges of 320-450 and 480-720 GeV. In Fig. 9, the NLO theoretical cross section for the gluon fusion production of a radion with decay constant Λ = 3 TeV is shown. Such a radion is excluded at 95% CL in the mass ranges of 260-280, 300-450, and 480-1120 GeV.

Summary
A search for a narrow-width resonance decaying into two Higgs bosons, each decaying into a bottom quark-antiquark pair, is presented. The search is performed using proton-proton collision data corresponding to an integrated luminosity of 35.9 fb −1 at √ s = 13 TeV recorded by the CMS detector at the LHC. No evidence for a signal is observed and upper limits at 95% confidence level on the production cross section for spin-0 and spin-2 resonances in the mass range from 260 to 1200 GeV are set. These cross-section limits are translated into an exclusion at 95% confidence level of a bulk KK-graviton (with κl = 35 and κ = 0.5M Pl ) in the mass ranges of 320-450 GeV and 480-720 GeV. The corresponding excluded mass ranges for a radion (with decay constant Λ = 3 TeV) are 260-280 GeV, 300-450 GeV and 480-1120 GeV. This analysis outperforms a similar search by CMS using 17.9 fb −1 collected at 8 TeV [15] and extends the sensitivity to the gluon fusion production of a radion with decay constant Λ = 3 TeV and to bulk graviton with κ set to 0.     [20] ATLAS Collaboration, "Search for Higgs boson pair production in the γγbb final state using pp collision data at √ s = 8 TeV from the ATLAS detector", Phys. Rev. Lett. 114 (2015) 081802, doi:10.1103/PhysRevLett.114.081802, arXiv:1406.5053.