Search for production of Higgs boson pairs in the four b quark final state using large-area jets in proton-proton collisions at $\sqrt{s} =$ 13 TeV

A search is presented for pair production of the standard model Higgs boson using data from proton-proton collisions at a centre-of-mass energy of 13 TeV, collected by the CMS experiment at the CERN LHC in 2016, and corresponding to an integrated luminosity of 35.9 fb$^{-1}$. The final state consists of two b quark-antiquark pairs. The search is conducted in the region of phase space where one pair is highly Lorentz-boosted and is reconstructed as a single large-area jet, and the other pair is resolved and is reconstructed using two b-tagged jets. The results are obtained by combining this analysis with another from CMS looking for events with two large jets. Limits are set on the product of the cross sections and branching fractions for narrow bulk gravitons and radions in warped extra-dimensional models having a mass in the range 750-3000 GeV. The resulting observed and expected upper limits on the non-resonant Higgs boson pair production cross section correspond to 179 and 114 times the standard model value, respectively, at 95% confidence level. The existence of anomalous Higgs boson couplings is also investigated and limits are set on the non-resonant Higgs boson pair production cross sections for representative coupling values.


Introduction
Models with a warped extra dimension (WED), as proposed by Randall and Sundrum [7], are among those BSM scenarios that predict the existence of resonances with large couplings to the SM Higgs boson, such as the spin-0 radion [9][10][11] and the spin-2 first Kaluza-Klein (KK) excitation of the graviton [12][13][14]. The WED models postulate an additional spatial dimension l compactified between two four-dimensional hypersurfaces known as the branes, with the region between, the bulk, warped by an exponential metric κl, where κ is the warp factor [15]. A value of κl∼35 fixes the mass hierarchy between the Planck scale M Pl and the electroweak scale [7]. One of the parameters of the model is κ/M Pl , where M Pl ≡ M Pl / √ 8π. The ultraviolet cutoff scale of the model Λ R ≡ √ 6e −κl M Pl [9] is another parameter, and is expected to be near the TeV scale.
In the absence of new resonances coupling to the Higgs boson, the gluon fusion Higgs boson pair production subprocess can still be enhanced by BSM contributions to the coupling parameters of the Higgs boson and the SM fields [16]. The SM production rate of HH through gluon fusion is determined by the Yukawa coupling of the Higgs boson to the top quark y SM t and the Higgs boson self-coupling λ SM HH = m 2 H /2v 2 . Here, m H = 125 GeV is the Higgs boson mass [17,18] and v = 246 GeV is the vacuum expectation value of the Higgs field. Deviations from the SM values of these two coupling parameters can be expressed as κ λ ≡ λ HH /λ SM HH and κ t ≡ y t /y SM t , respectively. Depending on the BSM scenario, other couplings not present in the SM may also exist and can be described by dimension-6 operators in the framework of an effective field theory by the Lagrangian [19]: The anomalous couplings and the corresponding parameters in this Lagrangian are: the contact interaction between a pair of Higgs bosons and a pair of top quarks (c 2 ), the interaction between the Higgs boson and the gluon (c g ), and the interaction between a pair of Higgs bosons and a pair of gluons (c 2g ). The couplings with CP-violation and the interactions of the Higgs boson with light SM and BSM particles are not considered. The Lagrangian models the effects of BSM scenarios with a scale that is beyond the direct LHC reach. This five-parameter space of BSM Higgs couplings has constraints from measurements of single Higgs boson production and other theoretical considerations [20,21].
Searches for HH production have been performed by the ATLAS [22][23][24][25][26] and CMS [27][28][29][30][31][32][33][34][35] Collaborations using the LHC pp collision data at √ s = 8 and 13 TeV. A search targeting the high m X range for a KK bulk graviton or a radion decaying to HH, in the bbbb final state, was published by the CMS Collaboration [36], in which two large-area jets are used to reconstruct the highly Lorentz-boosted Higgs bosons ("fully-merged" event topology). A similar search, focusing on a lower range of m X , was also performed by CMS [37], using events with four separate b quark jets. The configuration of a Higgs boson candidate as one large-area jet or as two separate smaller jets is dependent on the momentum of the Higgs boson [38].
In this paper, we improve upon the CMS search for high mass resonance (750 ≤ m X ≤ 3000 GeV) decaying to HH → bbbb [36] by using "semi-resolved" events, i.e. those containing exactly one highly Lorentz-boosted Higgs boson while the other Higgs boson is required to have a lower boost. The more boosted Higgs boson is reconstructed using a large-area jet and the other is reconstructed from two separate b quark jets. The inclusion of the semi-resolved events leads to a significant improvement in the search sensitivity for resonances with 750 ≤ m X ≤ 2000 GeV. With the addition of the semi-resolved events, a signal from the non-resonant production of HH is also accessible using boosted topologies, since such production typically results in an HH invariant mass that is lower than that of a postulated resonance signal. For full sensitivity, the results are obtained using a statistical combination of the semi-resolved events with the fully-merged events selected using the criteria in Ref. [36]. In addition to improving the search for X → HH, strong constraints are thus obtained for several regions in the H boson anomalous coupling parameter space, defined by Eq. (1).

The CMS detector and event reconstruction
The CMS detector with its coordinate system and the relevant kinematic variables is described in Ref. [39]. The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the field volume are silicon pixel and strip trackers, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. The tracker covers a pseudorapidity η range from −2.5 to 2.5 with the ECAL and the HCAL extending up to |η| = 3. Forward calorimeters in the region up to |η| = 5 provide good hermeticity to the detector. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid, covering a region of |η| < 2.4.
Events of interest are selected using a two-tiered trigger system [40]. The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage. Events used in this analysis are selected at the trigger level based on the presence of jets in the detector. The level-1 trigger algorithms reconstruct jets from energy deposits in the calorimeters. The particle-flow (PF) algorithm [41], aims to reconstruct and identify each individual particle in an event. The physics objects reconstructed include jets (clustered with a different algorithm), electrons, muons, photons, and also the missing-p T vector.
Multiple pp collisions may occur in the same or adjacent LHC bunch crossings (pileup) and contribute to the overall event activity in the detector. The reconstructed vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [42,43] with the tracks assigned tot he vertex as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets. The other interaction vertices are designated as pileup vertices.
The energy of each electron is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The energy of each muon is obtained from the curvature of the corresponding track. The energy of each charged hadron is determined from a combination of its momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers. Finally, the energy of each neutral hadron is obtained from the corresponding corrected ECAL and HCAL energies.
Particles reconstructed by the PF algorithm are clustered into jets with the anti-k T algorithm [42,43], using a distance parameter of 0.8 (AK8 jets) or 0.4 (AK4 jets). The jet transverse momentum is determined as the vector sum p T of all clustered particles. To mitigate the effect of pileup on the AK4 jet momentum, tracks identified as originating from pileup vertices are discarded in the clustering, and an offset correction [44,45] is applied for remaining contributions from neutral particles. Jet energy corrections are derived from simulation to bring the measured response of the jets to that of particle level jets on average. In situ measurements of the momentum balance in events containing either a pair of jets, or a Z boson or a photon recoiling against a jet, or several jets, are used to account for any residual differences in jet energy scale in data and simulation. Additional selection criteria are applied to each jet to remove jets potentially dominated by anomalous contributions from various subdetector components. After all calibrations, the jet p T is found from simulation to be within 5-10% of the true p T of the clustered particles, over the measured range [45,46].
For the AK8 jet mass measurement, the "pileup per particle identification" algorithm [47] (PUPPI) is applied to remove pileup effects from the jet. Particles from the PF algorithm, with their PUPPI weights, are clustered into AK8-PUPPI jets which are groomed [48] to remove soft and wide-angle radiation using the soft-drop algorithm [49,50], using the soft radiation fraction parameter z = 0.1 and the angular exponent parameter β = 0. Dedicated mass corrections [51], derived from simulation and data in a region enriched with tt events containing merged W → qq decays, are applied to the jet mass in order to remove residual dependence on the jet p T , and to match the jet mass scale and resolution observed in data. The AK8 jet soft-drop mass is assigned by matching the groomed AK8-PUPPI jet with the original jet using the criterion ∆R(AK8 jet, AK8-PUPPI jet) < 0.8, where ∆R ≡ √ (∆η) 2 + (∆φ) 2 , φ being the azimuthal angle in radians. Approximately 100% of the AK8 jets can be matched to an AK8 PUPPI jet.

Event simulation
The bulk graviton and radion signal events are simulated at leading order in the mass range 750-3000 GeV with a width of 1 MeV (much smaller than experimental resolution), using the MADGRAPH5 aMC@NLO 2.3.3 [52] event generator. The NNPDF3.0 leading order parton distribution function (PDF) set [53], taken from LHAPDF6 PDF set [54][55][56][57], with the fourflavour scheme, is used. The showering and hadronization of partons are simulated with PYTHIA 8.212 [58].
The HERWIG++ 2.7.1 [59] generator is used as an alternative model, to evaluate the systematic uncertainty associated with the parton shower and hadronization. The tune CUETP8M1-NNPDF2.3LO [60] is used for PYTHIA 8, while the EE5C tune [61] is used for HERWIG++.
Non-resonant HH signals were generated using the effective field theory approach defined in Refs. [4,62] and is described by the five parameters given in Eq. 1: κ λ , κ t , c 2 , c g , and c 2g . The final state kinematic distributions of the HH pairs depend upon the values of these five parameters. A statistical approach was developed to identify twelve regions of the parameter space, referred to as clusters, with distinct kinematic observables of the HH system. In particular, models in the same cluster have similar distributions of the di-Higgs boson invariant mass m HH , the transverse momentum of the di-Higgs boson system, and the modulus of the cosine of the polar angle of one Higgs boson with respect to the beam axis, while the distributions of these variables are unique when comparing models from different clusters [63]. For each cluster, a set of representative values of the five parameters is chosen, referred to as the "shape benchmarks". Events are simulated for each of these shape benchmarks, as well as for the SM values of these couplings, and the case where the Higgs boson self-coupling vanishes, i.e. κ λ = 0. The values of these benchmark coupling parameters are given in Table 1. Table 1: Parameter values of the couplings corresponding to the twelve shape benchmarks, the SM prediction, and the case with vanishing Higgs boson self-interaction, κ λ = 0. The dominant background consists of events comprised uniquely of jets (multijet events) arising from the SM quantum chromodynamics (QCD) interaction, and is modelled entirely from data. The remaining background, consisting mostly of tt+jets events, is less than 10% of the total background, is modelled using POWHEG 2.0 [64][65][66] and interfaced to PYTHIA 8. The CUETP8M2T4 tune [67,68] is used for generating the tt+jets events. The tt+jets background rate is estimated using a next-to-next-to-leading order cross section of 832 +46 −52 pb [69], corresponding to the top quark mass of 172.5 GeV. A sample of multijet events from QCD interactions, simulated at leading order using MADGRAPH5 aMC@NLO and PYTHIA 8, is used to develop and validate the background estimation techniques, prior to being applied to the data.
All generated samples were processed through a GEANT4-based [70,71] simulation of the CMS detector. The effect of pileup, averaging to 23 at the LHC beam conditions in 2016, is included in the simulations, and the samples are reweighted to match the distribution of the number of pp interactions observed in the data, assuming a total inelastic pp collision cross section of 69.2 mb [72].

Event selection
Five different HLT triggers were used to collect the semi-resolved events used in this analysis. An event is selected if the scalar sum of the p T of all AK4 jets in the event (H T ) is greater than 800 or 900 GeV, depending on the LHC beam instantaneous luminosity. Events with H T ≥ 650 GeV, and a pair of jets with invariant mass above 900 GeV and a pseudorapidity separation |∆η| < 1.5 are also selected. A third HLT trigger accepts events if the scalar sum of the p T of all AK8 jets is greater than 650 or 700 GeV and the "trimmed mass" of an AK8 jet is above 50 GeV. The jet trimmed mass is obtained after removing remnants of soft radiation with the jet trimming technique [73], using a subjet size parameter of 0.3 and a subjet-to-AK8 jet p T fraction of 0.1. Should an event contain an AK8 jet with p T > 360 GeV and a trimmed mass greater than 30 GeV, it is selected by the fourth HLT trigger. Events containing two AK8 jets having p T > 280 and 200 GeV, with at least one having trimmed mass greater than 30 GeV together with an AK4 jet passing a loose b-tagging criterion, pass the fifth HLT trigger.
Jets in events collected using the logical OR of the above HLT triggers are required to have |η| < 2.4, and p T > 30 GeV for AK4 jets and p T > 300 GeV for AK8 jets. One AK8 jet is used to identify a boosted and spatially merged H → bb decay (H jets) while two AK4 jets are used to reconstruct a spatially resolved H → bb decay.
The first H-tagging criterion requires an AK8 jet to have a soft-drop mass m J between 105 and 135 GeV, consistent with the measured mass of the Higgs boson m H = 125 GeV. This selection corresponds to an efficiency of about 60-70% for a resonant signal mass m X in the range 750-3000 GeV. The soft-drop jet mass interval was chosen to include a large fraction of the boosted H → bb signal, while avoiding overlaps with CMS analyses searching for bulk gravitons and radions decaying to boosted W and Z bosons [74]. The "N-subjettiness" algorithm [75] is used on the AK8-PUPPI jet constituents, to compute the variables τ N , which quantify the degree to which a jet contains N subjets. A selection on the ratio τ 21 ≡ τ 2 /τ 1 < 0.55 is required for all AK8 jets to be H tagged, which has a jet p T -dependent efficiency of 50-70%. The selection criterion on τ 21 was optimized for signal sensitivity over the range of m X values explored.
A jet flavour requirement using a "double-b tagger" algorithm [76] is applied to the AK8 jet as the final H-tagging requirement. The double-b tagger is a multivariate discriminator with an output between -1 and 1, a higher value indicating a greater probability for the jet to contain a bb pair. The double-b tagger exploits the presence of two hadronized b quarks inside the boosted H → bb decay, and uses variables related to b hadron lifetime and mass to distinguish H jets against a background of jets of other flavours. The double-b tagger algorithm also exploits the fact that the directions of the b hadron momentum are strongly correlated with the axes used to calculate the N-subjettiness variables τ N . An H jet candidate should have a double-b tagger discriminator greater than 0.8, which corresponds to an efficiency of 30% and a misidentification rate of about 1%, as measured in a sample of multijet events. The efficiency of the double-b tagger for simulated jets is corrected to match that in the data, based on efficiency measurements using jets containing pairs of muons, thereby yielding samples enriched in jets from gluons splitting to bb pairs. These efficiency corrections are in between 0.92 and 1.02, for jets in the selected p T range.
To find a Higgs boson decay into two resolved b quark jets, all AK4 jets in each event are examined for their b tag value using "DeepCSV" algorithm, which is a deep neural network, trained using information from tracks and secondary vertices associated to the jets [76]. The DeepCSV discriminator gives the probability of a jet to have originated from the hadronization of a bottom quark. A selection on the DeepCSV discriminator of AK4 jets is made, corresponding to a 1% mistag rate for light flavoured jets. The corresponding b-tagging efficiency is about 70% for b quark jets in the p T range 80-150 GeV, and decreases to about 50% for p T ∼ 1000 GeV. The b tagging efficiency in the simulations is corrected to match the one in the data, using measurements of the b tagging algorithm performance in a sample of muon-tagged jets and b jets from tt+jets events, where the correction factor ranges from approximately 0.95 to 1.1.
To identify events with a resolved H → bb decay, all pairs of b-tagged AK4 jet are examined, to find events with at least one pair where each AK4 jet is at least ∆R > 0.8 away from the leading-p T AK8 jet and within ∆R < 1.5 of each other. If several such pairs are found, the pair of jets, j 1 and j 2 , that has the highest sum of the AK4 jet DeepCSV discriminator values is selected. The leading-p T AK8 jet is then identified as the boosted H candidate, and the pair of AK4 jets is identified as the resolved H candidate. If no pairs are found, this process is repeated with the subleading-p T AK8 jet. If a pair of AK4 jets is identified, then the subleading-p T AK8 jet is identified as the boosted H candidate, and the pair of AK4 jets is identified as the resolved H candidate. If no pairs are found once again, the event is rejected. The invariant mass of j 1 and j 2 , m jj (j 1 , j 2 ), is required to be within 90-140 GeV, forming the resolved H → bb candidate.
The tt+jets background is reduced by reconstructing a t → bqq system in events with three or more AK4 jets, combining j 1 and j 2 with the nearest AK4 jet j 3 . For the tt+jets background, the trijet invariant mass m jjj (j 1 , j 2 , j 3 ) peaks around the top quark mass of 172 GeV. Hence, m jjj (j 1 , j 2 , j 3 ) is required be greater than 200 GeV, namely above the top quark mass. Events containing leptons (electrons or muons) with p T > 20 GeV and |η| < 2.4, and only a small amount of energy in an area around the lepton direction compared to the lepton p T , are rejected, to further suppress tt+jets and other backgrounds.
A resonant HH signal results in a small pseudorapidity separation between the two Higgs bosons, while the candidates from the multijet background typically have a larger pseudorapidity separation. Events are therefore categorized according to the pseudorapidity difference between the H jet and the resolved H → bb candidate. These two categories are defined by |∆η(H jet, resolved H → bb)| within the interval 0.0-1.0 or within the interval 1.0-2.0.
To search for resonant and a non-resonant HH signals, the invariant mass distribution of the boosted and resolved Higgs boson candidate system (m Jjj ) is examined for an excess of events over the estimated background. The "reduced di-Higgs invariant mass" is defined as The quantity m Jjj,red is used rather than m Jjj since by subtracting the masses of the reconstructed H candidates and adding back the exact Higgs boson mass m H , fluctuations from the jet mass resolution are corrected, leading to 8-10% improvement in the HH mass resolution. After the full selection, the multijet background is about 90% of the total background, with the remaining background being tt+jets events.
With the above event selection, the trigger criteria reach an efficiency of greater than 99% for events with m Jjj,red ≥ 1100 GeV. For lower values of invariant mass (between 750 and 1100 GeV), the trigger efficiency is between 80 and 99% for 0 ≤ |∆η(H jet, resolved H → bb)| < 1 and between 60 and 99% for 1 ≤ |∆η(H jet, resolved H → bb)| ≤ 2. The trigger efficiency for the data is estimated from a multijets sample collected with a control trigger requiring a single AK4 jet with p T > 260 GeV. The trigger efficiency for the simulated samples is corrected using a scale factor to match the observed efficiency for the data. This scale factor depends mildly on |∆η(H jet, resolved H → bb)|, and is hence applied as a function of this variable.
The AK8 jet soft-drop mass distribution, the N-subjettiness ratio τ 21 distribution, and the doubleb tagger discriminator distribution for the backgrounds and simulated signals are shown in Fig. 1. The DeepCSV discriminator distributions for the two AK4 jets, the dijet invariant mass distribution, and the trijet invariant mass distribution for the backgrounds and simulated sam-  ples are shown in Fig. 2. The selection criteria for the above plots is as follows: AK8 jets with p T > 300 GeV, AK4 jets with p T > 30 GeV, AK8 and AK4 jets with |η| < 2.4, AK8 jet soft-drop mass > 40 GeV, AK4 jets DeepCSV discriminator > 0.2219, ∆R < 1.5 separation between the AK4 jets, and ∆R > 0.8 separation between the AK8 jet and each AK4 jet.
The semi-resolved event selection is summarized in Table 2, where in addition to these criteria, the events that are used by the fully-merged analysis of Ref.
[36] are removed, as detailed at the end of this section. The event selection efficiencies for bulk gravitons and radions are given in Fig. 3, for different assumed masses in the range 750-2000 GeV. The selection efficiencies for the non-resonant signals are between 0.01-2%.
In view of the statistical combination of the semi-resolved and the fully-merged analyses, we briefly describe the search in the fully-merged topology [36]. The analysis in the fully-merged regime uses the same trigger selection and the same selection for the H jet identification, except for a different requirement on the double-b tagger. These events contain two H jets J 1 and J 2 instead of one. The fully-merged events are classified according to the values of the double-b tagger discriminators of the two H jets, with both J 1 and J 2 required to pass a loose double-b tagger discriminator value of > 0.3. Events are then categorized into those with both J 1 and J 2 passing a tighter double-b tagger discriminator requirement of > 0.8, and the rest. The pseudorapidity separation between J 1 and J 2 is required to be |∆η(J 1 , J 2 )| < 1.3. The reduced   A Higgs boson candidate which passes the boosted AK8 jet selection can also pass the selection for two resolved AK4 jets. In particular, signal samples with higher mass that pass the semi-resolved selection often pass the fully-merged selection because both Higgs candidates are merged, but one candidate still passes the selection for a resolved jet as well. For each signal, the final semi-resolved selection includes anywhere from 23-53% events that are used by the fully-merged analysis, whether in the signal region or to estimate the QCD multijets background. These events are then removed from the semi-resolved analysis to allow for a combination with the fully-merged analysis.

Background estimation
The multijet background estimation technique for the semi-resolved analysis is the same as that for the fully-merged analysis [36]. A set of signal-free control regions is defined by changing the criteria on the soft-drop mass and the double-b tagger discriminator of the selected AK8 jet from those used for the H tagging. The selection criteria applied to the AK4 jets forming  Figure 3: The signal selection efficiencies for the radion and the bulk graviton, for different masses. The events are required to pass the selections given in Table 2 as well as to fail the selections of the fully-merged analysis of Ref. [36].
the resolved H → bb are the same as those used for the signal regions. If the soft-drop mass is within 60 GeV above or below the H jet mass window of 105-135 GeV, these regions are referred to as the mass sideband regions. These sidebands are separated into regions that pass or fail the double-b tagger tagging requirement.
We define the pass-fail ratio R p/f as the ratio of events for which the AK8 jet passes and fails the double-b tagger tagging requirement. The R p/f is measured in the soft-drop mass sidebands as a function of soft-drop mass. These values are fit to a quadratic function of the H jet mass to calculate the R p/f in the signal region. The antitag region, defined with the same criteria as the signal region, but with the AK8 jet failing the double-b tagger requirement, is then scaled by the R p/f value to estimate the multijets background in the signal region. This is done in bins of soft-drop mass to predict the entire background shape. The dependence of R p/f on m Jjj,red was found to be negligible, within the measurement uncertainties. Both the shape of the background m Jjj,red distribution and its total yield in the signal region is obtained using this method.
Prior to estimating the background, the tt+jets contributions derived from Monte Carlo simulation are subtracted from all sideband and signal regions in the data, and then added back in . The fitted function is interpolated to obtain R p/f in the signal region. Lower: The reduced mass distribution m Jjj,red in the data (black markers) with the estimated background represented as the black histogram. The tt+jets contribution from simulation is represented in green. The rest of the background is multijets, calculated by applying the R p/f to the antitag region. The uncertainty in the total background, before fitting the background model to the data, is depicted using the shaded region. The signal distributions for a bulk graviton with a mass of 800 GeV (blue) and the non-resonant benchmark 2 model (red) are also shown for assumed values of the products of the production cross sections for HH and the branching fraction to 4b, σB. For the upper and lower figures, the pseudorapidity intervals are 0 ≤ |∆η| < 1 and 1 ≤ |∆η| ≤ 2, respectively. once the multijet background calculation is completed, to estimate the contribution of tt+jets to the total background. The fractions of signal events in the sideband regions were found to be negligible as compared with the total numbers of events. Figure 4 (left) shows the quadratic fit in the AK8 jet soft-drop mass sidebands of the pass-fail ratio R p/f as a function of AK8 jet soft-drop mass, as obtained in the data and in the predicted background shape in the signal region, where overlap with the merged analysis in the signal, sideband, and antitag regions is removed. The functional form was chosen after performing a Fisher F-test [77], which established that, among polynomials, a quadratic form is necessary and sufficient. Other functional forms were tested and the fit results were found to be consistent with that using the quadratic function. The resulting background distributions are compared with the observed data, as shown in Fig. 4 (right).

Systematic uncertainties
The following sources of systematic uncertainty affect the expected signal and background event yields. None of these lead to a significant change in the signal shape. A complete list of systematic uncertainties is given in Table 3. The trigger response modelling uncertainties are particularly important for m Jjj,red < 1100 GeV, where the trigger efficiency drops below 99%. The trigger efficiency data-to-simulations scale factor has an uncertainty between 1 and 15%, attributable to the control trigger inefficiency and the sample size used.
The impact of the jet energy scale and resolution uncertainties [51] on the signal yields was estimated to be 1-3%, depending on the signal mass. The jet mass scale and resolution, as well as the τ 21 selection efficiency data-to-simulation scale factors were measured using a sample of boosted W → qq jets in semileptonic tt events. The jet mass scale and resolution has a 2% effect on the signal yields because of a change in the mean of the H jet mass distribution. A correction factor is applied to account for the difference in the jet shower profile of W → qq and H → bb decays, by comparing the ratio of the efficiency of H and Wjets using the PYTHIA 8 and HERWIG++ shower generators. This uncertainty, the H tagging correction factor, is in the range 5-20%, depending on the resonance mass m X . The τ 21 selection efficiency uncertainty depends on how many τ 21 tags are used, two for the fully-merged (26-30% uncertainty) and one for the semi-resolved analysis (13-15% uncertainty). This includes an additional uncertainty in the τ 21 scale factor, determined using simulations, for jets with p T higher than those in the tt events used for the evaluation of this systematic.
Scale factors are used to correct the signal events yields so their double-b tagger and DeepCSV discriminator efficiencies are the same as for data. The double-b tagger and the DeepCSV discriminator scale factors are taken to be 100% correlated. The associated uncertainty is 2-9% [76], depending on the double-b tagger and requirement threshold and jet p T , and is propagated to the total uncertainty in the signal yield.
The impact of the theoretical scale uncertainties and PDF uncertainties, the latter derived using the PDF4LHC procedure [57] and the NNPDF3.0 PDF sets, is estimated to be 0.1-3%. These uncertainties affect the product of the signal acceptance and the selection efficiency. The scale and the PDF uncertainties have negligible impact on the signal m Jjj,red distributions. Additional systematic uncertainties associated with the pileup modelling (1-2%, based on a 4.6% variation on the pp total inelastic cross section) and with the integrated luminosity determination (2.5%) [78], are applied to the signal yield.
The systematic uncertainty on the trijet invariant mass cut was calculated by comparing the cut efficiency for Pythia and Herwig bulk graviton samples, and is equivalent to a 0.5% systematic.
The systematic uncertainty applied to the signal is also applied to the tt+jets background in the semi-resolved analysis, as appropriate. The total uncertainty in the tt+jets background is 11-15%, of which 6% derives from the uncertainty in the tt+jets cross section.
The main source of uncertainty for the multijet background is due to the statistical uncertainty in the fit to the R p/f ratio performed in the H jet mass sidebands. This uncertainty, amounting to 2-10%, is fully correlated between all m Jjj,red bins. Additional statistical uncertainties on the background shape and yield in the signal region result from the finite statistics of the multijets samples in the antitag region and are evaluated using the Barlow-Beeston Lite method [79,80]. These uncertainty are small as compared with the uncertainty on the R p/f ratio, and are uncorrelated from bin to bin.

Results
This analysis extends the search for a resonance X decaying to HH → bbbb with two boosted H jets [36] to cover the semi-resolved topology involving one boosted H jet and one resolved H → bb decay reconstructed using two b jets. An HH signal would appear as an excess of events over estimated background in the m Jjj,red spectra of the different signal event categories, as discussed in Section 5.
The binned m Jjj,red distributions of the signal and the background are fitted to the data, as shown in Fig. 4 (right), and examined for an excess of events above the predicted background. The data were found to be consistent with the expected background predictions. Upper limits on the product of the signal cross sections and branching fractions are obtained using the profile likelihood as a test statistic [81]. The systematic uncertainties are treated as nuisance parameters and are profiled in the minimization of the negative of the logarithm of the profile likelihood ratio and the distributions of the likelihood ratio are calculated using the asymptotic approximation [82] of the procedure reported in Refs. [83,84]. Upper limits at 95% confidence level are set on the product of the production cross section and the branching fractions σ(pp → X)B(X → bbbb).
Results are obtained using a statistical combination of the semi-resolved and fully-merged event categories for the bulk graviton having a mass between 750-2000 GeV, and a radion with a mass between 750-1600 GeV. Above these mass ranges, the inclusion of the semi-resolved events does not appreciably improve the search sensitivity. The limits on σ(pp → X)B(X →    bbbb) are shown in Fig. 5, and tabulated in Tables 4 and 5 for the bulk graviton and the radion, respectively. The limits for m X > 2000 GeV for the bulk graviton, and m X > 1600 GeV for the radion are those from the fully-merged analysis of Ref. [36].
For the interpretation of the results, this paper uses the scenario of Ref. [85] to describe a KK graviton, where the propagation of SM fields is allowed in the bulk, and follows the characteristics of the SM gauge group, with the right-handed top quark localized near the TeV brane.
The radion is an additional element of WED models that is needed to stabilize the size of the extra dimension l. The theoretical cross sections for σ(pp → X)B(X → bbbb) are calculated using κ/M Pl = 0.5 for the bulk gravitons and Λ R = 3 TeV for the radions, of different masses.
For these values of κ/M Pl and Λ R , the branching fractions B(X → bbbb) are 10 and 23%, for the graviton and the radion, respectively [86]. As shown in Fig. 5 (right), a radion having a mass between 1000 and 1500 GeV is excluded at 95% confidence level for Λ R = 3 TeV.
The improvement in the upper limits on σ(pp → X)B(X → bbbb) due to the inclusion of the semi-resolved event category between 18 and 7%, for a bulk graviton in the mass range 750-2000 GeV. A much larger improvement-between 55 and 8%-is seen for a radion in the mass range 750-1600 GeV. This can be attributed to the two pseudorapidity intervals, |∆η| < 1 and 1 ≤ |∆η| ≤ 2, utilized in the semi-resolved event selection, with the lower pseudorapidity interval having a better signal to background ratio for a spin-0 radion, because of the angular distribution of its decay products. In addition, both the fully-merged and the semi-resolved analyses look for non-resonant HH production. The observed and expected upper limits are presented in Table 6 for the semiresolved and the fully-merged signal categories combined, also depicted in Fig. 6. The observed and expected limits are respectively, 179 and 114 times the product of the SM cross sections and branching fractions. The new limits are better by about a factor of three for benchmark 2 and a factor of two for benchmark 5, with significant improvements for benchmarks 8, 9, and 11, compared to existing measurements [33,35]. The increased sensitivity of this analysis to certain benchmarks is due to the higher level of destructive interference among the HH production processes close to the kinematic threshold, which leads to a corresponding shift of the HH mass spectrum towards higher values. This leads to, on average, a higher p T of the Higgs bosons, and hence in the sensitivity of this analysis, which identifies Higgs bosons using boosted techniques.

Summary
A search is presented for the pair production of standard model Higgs bosons (HH), both decaying to a bottom quark-antiquark pair (bb), using data from proton-proton collisions at  The results of the search are compared with predictions for the resonant production of a narrow Kaluza-Klein bulk graviton and a narrow radion in warped extradimensional models. The search is also sensitive to several beyond standard model non-resonant HH production scenarios. Such cases may arise either when an off-shell massive resonance produced in proton-proton collisions decays to HH, or through beyond standard model effects in the Higgs boson coupling parameters. The results are interpreted in terms of upper limits on the product of the cross section for the respective signal processes and the branching fraction to HH → bbbb, at 95% confidence level.
The upper limits range from 43.9 to 1.4 fb for the bulk graviton and from 67 to 1.6 fb for the radion for the mass range 750-3000 GeV. Depending on the mass of the resonance, these limits improve upon the results of Ref.
[36] by up to 18% for the bulk graviton and up to 55% for the radion.
The non-resonant production of Higgs boson pairs is modelled using an effective Lagrangian with five coupling parameters. The upper limit corresponding to the standard model values of the coupling parameters is placed at 1980 fb, which is 179 times the prediction. In addition, upper limits in the range of 4520 to 36.7 fb are set on twelve shape benchmarks, i.e. representative sets of the five coupling parameters [63]. These are the first limits on non-resonant Higgs boson pair-production signals using boosted topologies, and are the most stringent limits to date for the shape benchmarks 2, 5, 8, 9, and 11.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses.   [21] ATLAS Collaboration, "Measurements of the Higgs boson production and decay rates and coupling strengths using pp collision data at √ s = 7 and 8 TeV in the ATLAS experiment", Eur. Phys. J. C 76 (2016) 6, doi:10.1140/epjc/s10052-015-3769-y, arXiv:1507.04548. [24] ATLAS Collaboration, "Searches for Higgs boson pair production in the HH → bbττ, γγWW * , γγbb, bbbb channels with the ATLAS detector", Phys. Rev. D 92 (2015) 092004, doi:10.1103/PhysRevD.92.092004, arXiv:1509.04670.
[26] ATLAS Collaboration, "Search for pair production of Higgs bosons in the bbbb final state using proton-proton collisions at √ s = 13 TeV with the ATLAS detector", (2018). arXiv:1804.06174.