Search for resonant pair production of Higgs bosons in the bbZZ channel in proton-proton collisions at $\sqrt{s}=$ 13 TeV

A search for the production of a narrow-width resonance decaying into a pair of Higgs bosons decaying into the bbZZ channel is presented. The analysis is based on data collected with the CMS detector during 2016, in proton-proton collisions at the LHC, corresponding to an integrated luminosity of 35.9 fb$^{-1}$. The final states considered are the ones where one of the Z bosons decays into a pair of muons or electrons, and the other Z boson decays either to a pair of quarks or a pair of neutrinos. Upper limits at 95% confidence level are placed on the production of narrow-width spin-0 or spin-2 particles decaying to a pair of Higgs bosons, in models with and without an extended Higgs sector. For a resonance mass range between 260 and 1000 GeV, limits on the production cross section times branching fraction of a spin-0 and spin-2 resonance range from 0.1 to 5.0 pb and 0.1 to 3.6 pb, respectively. These results set limits in parameter space in bulk Randall-Sundrum radion, Kaluza-Klein excitation of the graviton, and N2HDM models. For specific choices of parameters the N2HDM can be excluded in a mass range between 360 and 620 GeV for a resonance decaying to two Higgs bosons. This is the first search for Higgs boson resonant pair production in the bbZZ channel.


Introduction
The discovery of the Higgs boson (h) in 2012 [1][2][3] has led to a detailed program of studies of the Higgs field couplings to the elementary particles of the standard model (SM) of particle physics: leptons, quarks, and gauge bosons. To fully understand the form of the Higgs field potential, which is a key element in the formulation of the SM, it is important to also study the self-interaction of the Higgs boson. The self-interaction can be investigated through measurements of the production of a pair of Higgs bosons (hh). In the SM, hh production is a rare, nonresonant process, with a small production rate [4] that will require the future data sets of the high-luminosity LHC to be observed [4]. Hence, an early observation of hh production, a resonant production in particular, would be a spectacular signature of physics beyond the standard model (BSM). The production of gravitons, radions, or stoponium [5][6][7], for example, could lead to s-channel hh production via narrow-width resonances. The breadth of the Higgs boson decay channels provides a unique opportunity to test the self-consistency of an hh signal with the SM or models with extended electroweak sectors, such as two-Higgs doublet models (2HDM) [8,9] or extensions of the minimal supersymmetric standard model [10][11][12].
This paper reports a search for resonant pp → X → HH production in the HH → bbZZ decay channel, where X is a narrow-width resonance of spin-0 or spin-2, and H can represent either h or an additional Higgs boson from an extended electroweak sector. The search uses proton-proton (pp) collision data at √ s = 13 TeV, recorded with the CMS detector at the LHC in 2016, and corresponding to an integrated luminosity of 35.9 fb −1 . It covers a range of resonance masses between 260 and 1000 GeV. The final state consists of two b jets from one Higgs boson decay and two distinct Z boson decay signatures from the other H → ZZ decay: two same-flavor, opposite-sign (OS) leptons from a decay of one of the Z bosons, and either two jets of any flavor (the bb jj channel) or significant missing transverse momentum (the bb νν channel) from the decay of the second Z boson to neutrinos. In both cases, the selected charged leptons are either electrons or muons. In the SM, the branching fractions of these signatures represent 0.43 (0.12)% of the full hh decay through the bbZZ intermediate state in the bb jj (bb νν) channel. The challenging aspect of the search in the bb jj channel is the ability to discriminate the signal containing two b jets and two additional jets from multijet background events. For a search in the bb νν channel, the challenge resides in discriminating the signal against top quark anti-top quark (tt) events and instrumental background sources of large missing transverse momentum arising from the mismeasurement of the energies of jets in the final state. The two channels are kept independent by applying orthogonal selections on the missing transverse momentum of the event. Signal yields are calculated for each individual channel and are then combined. Having multiple decay channels with complementary background compositions and sensitivities over a large resonance mass (m X ) range makes this combination of the bb νν and bb jj channels highly efficient for covering the bbZZ final state. This is the first search for Higgs boson resonant pair production in the bbZZ channel.
Previous searches for resonant hh production have been performed by the CMS and ATLAS Collaborations in the bbbb [13,14], bbττ [15,16], bbγγ [17], and bb ν ν [16,18] channels. While coverage of as many hh decay channels as possible remains necessary to understand the exact nature of the Higgs boson self-coupling and the electroweak symmetry breaking mechanism, a bbZZ search is particularly interesting in models with extended electroweak sectors, where the phenomenology of additional Higgs bosons can lead to significantly enhanced bbZZ production, while suppressing the BSM production of bbbb, bbττ, or bbγγ final states.

Benchmark models
As in the previous searches, a class of narrow width resonance models arising from the Randall-Sundrum (RS) model [19] in warped extra dimensions [20][21][22][23] are considered. This scenario introduces one small spatial extra dimension with a nonfactorizable geometry, where the SM particles are not allowed to propagate along that extra dimension, and is referred to in this search as RS1. The resonant particle can be a radion (spin-0) or the first Kaluza-Klein (KK) excitation of a graviton (spin-2). The production cross section of the radion is proportional to 1/λ 2 R where λ R is the interaction scale parameter of the theory. In this analysis, we consider the cases where λ R = 1 TeV with kL = 35, where k is the constant in the warp factor (e −kL ) appearing in the space-time metric of the theory and L is the size of the extra dimension. The free parameter of the model for the graviton case isk = k/M Pl , where M Pl is the reduced Planck scale, and we considerk = 0.1 in this analysis [24]. We further scan the model parameter space in the λ R andk parameters for their respective models. Production at hadron colliders is expected to be dominated by gluon-gluon fusion, and we assume that the radion or graviton is produced exclusively via this process. Due to the small branching fraction of hh → bbZZ and the high multiplicities of the final states, the analyses presented in this paper are less sensitive to these models compared to the previous searches. As noted in Section 1, however, certain models with extended electroweak sectors can produce significantly enhanced bbZZ production, while suppressing final states with Higgs boson decays to fermions and scalar bosons.
Such an enhancement can be produced for example in the next-to-minimal 2HDM (N2HDM) extended Higgs sector [25,26], where an additional real singlet is introduced in addition to the usual two doublet Higgs bosons of the 2HDM. This analysis is further interpreted in this context. The so-called broken phase is considered, wherein both the Higgs doublets and the singlet acquire vacuum expectation values (vev) [26]. Mixing between these states produces 3 CP-even Higgs bosons H 1 , H 2 , and H 3 , with masses that are free parameters of the model. This search considers the nearly mass-degenerate case where the masses of the two bosons H 1 and H 2 are constrained to the experimental measurements of the h mass, which would be indistinguishable from h production with current LHC data sets [10,27,28], but may give rise to manifestly non-SM-like rates in the case of hh production. In what is commonly referred to as Higgs boson cascade decays, H 3 can decay to any combination of bosons H 1 and H 2 , which then both can have different decay branching fractions compared to the SM Higgs boson. The model spectrum depends on the ratio of the vevs of the two Higgs doublets tan β, low values of which enhance H 3 production; the vev of the singlet, which affects the decay branching fractions of H 3 to H 1 and H 2 ; and three mixing angles, which determine the decay branching fractions of H 1 and H 2 [26]. The model spectra described below are determined using N2HDECAY [29], and are chosen to enhance production of the bbZZ final state while roughly respecting current LHC measurements of the SM h branching fractions [4]. The gluon-gluon fusion production cross sections of H 3 are determined from the BSM Higgs boson predictions of the LHC Higgs Cross Section Working Group [4]. These cross sections assume SM decay branching fractions of the Higgs boson, and changing these branching fractions affects the production cross section. The cross sections are corrected at leading order (LO) by the ratio of the relative partial width of H 3 in the decay to two gluons compared to the BSM Higgs boson prediction. Enhanced (reduced) coupling of H 3 to gluons will enhance (reduce) the production cross section of H 3 . The mass of the Higgs bosons H 1 and H 2 are set to 125 GeV, and the mass of H 3 is generated in the range 260 ≤ m X ≤ 1000 GeV. Two benchmark points are chosen for this analysis, corresponding to tan β = 0.5 and 2.0. In both cases, the scalar vev is set to 45 GeV, and the mixing angles α 1 , α 2 , α 3 are set to 0.76, 0.48, and 1.00, respectively. For tan β = 0. 5 . This represents a 5% increase in the branching fraction to bbZZ compared to SM hh decays. The correction factor based on the relative partial width of H 3 to two gluons is around 0.7. These corrections and branching fractions produce significant differences in the production rates of the bbZZ system compared to hh production both in the SM and through resonant production of radions or gravitons.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity coverage provided by the barrel and endcap detectors, where pseudorapidity is defined as η = − ln[tan(θ/2)], and θ is the polar angle. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. CMS uses a two-level trigger system [30]. The first level of the CMS trigger system, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select the most interesting events. The high-level trigger (HLT) processor farm further decreases the event rate from around 100 kHz to a rate of around 1 kHz, before data storage. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [31].

Event simulation
The signal samples of RS1 spin-0 radion and RS1 KK spin-2 graviton narrow resonances decaying to a pair of Higgs bosons (X → hh) are generated at LO using MADGRAPH5 aMC@NLO. The h mass is set to 125 GeV, and the m X is generated in the range of 260-1000 GeV. In the bb νν channel the final state can be produced via either the bbZZ or bbW + W − intermediate states.
The main background processes to production of a pair of Higgs bosons in the bbZZ → bb jj or bb νν final states are Z/γ * +jets and tt processes. Less significant backgrounds arise from single top quark, W+jets, diboson+jets, SM Higgs boson production, and quantum chromodynamics (QCD) multijet production. Signal and background processes are modeled with simulations, with the exception of the QCD multijet background that is estimated using data control regions.
In the analysis using the bb jj channel, the Z/γ * +jets and W+jets processes are generated with MADGRAPH5 aMC@NLO2.4.2 [32] at next-to-leading order (NLO). In this case, the generator uses the FXFX jet merging scheme [33]. The analysis of the bb νν channel uses samples of Z/γ * +jets generated with MADGRAPH5 aMC@NLO at LO, with the MLM matching scheme [34], and reweighted to account for higher order QCD and electroweak effects [35].
The simulated samples are normalized to their best-known highest-order-QCD cross sections, either evaluated at NLO with MCFM [42] (diboson+jets) or at next-to-next-to-leading order with FEWZ 3.1 [43] (single top quark, W+jets, SM Higgs boson), with the exception of tt and Z/γ * +jets processes, which are normalized using data.
The simulated samples are interfaced with PYTHIA 8.212 [44] for parton showering and hadronization. The PYTHIA generator uses the CUETP8M1 underlying event tune [45]. The NNPDF3.0 NLO and LO parton distribution functions (PDFs) [46] are used for the various processes, with the precision matching that in the matrix element calculations.
For all the simulated samples used in this analysis, a simulation of CMS detector response based on GEANT4 [47] is applied. The presence of additional interactions in the same bunch crossing (pileup, or PU), both in-time and out-of-time with respect to the primary interaction, is simulated and corrected to agree with a multiplicity corresponding to the distribution measured in data.

Event reconstruction
Events are selected using triggers that require two muons with transverse momentum p T > 17 (8) GeV or two electrons with p T > 23 (12) GeV for the leading (sub-leading) lepton.
The particle-flow (PF) algorithm [48], which combines information from various elements of the CMS detector, is used to reconstruct and identify final-state particles, such as photons, electrons, muons, and charged and neutral hadrons, as individual PF objects. Combinations of PF objects are then used to reconstruct higher-level objects such as jets and missing transverse momentum.
Jets are reconstructed from the PF objects, using the anti-k T [49,50] algorithm with a distance parameter of R = 0.4. In order to reduce instrumental backgrounds and the contamination from PU, selected jets are required to satisfy loose identification criteria [51] based on the multiplicities and energy fractions carried by charged and neutral hadrons. The energy of reconstructed jets is calibrated using p T -and η-dependent corrections to account for nonuniformity and nonlinearity effects of the ECAL and HCAL energy response to neutral hadrons, for the presence of extra particles from PU, for the thresholds used in jet constituent selection, reconstruction inefficiencies, and possible biases introduced by the clustering algorithm. These jet energy corrections are extracted from the measurement of the momentum balance in dijet, photon + jet, Z/γ * +jets, and multijet events [52]. A residual η-and p T -dependent calibration is applied to correct for the small differences between data and simulated jets. The jets that are candidates to be from the decay of one of the Higgs bosons and of one of the Z bosons are required to have p T > 20 GeV. Furthermore, jets are required to have a spatial separation of ∆R > 0.3 from lepton candidates.
Jets originating from b quarks are identified with the combined multivariate analysis (cMVA) algorithm [53]. A jet is tagged as a b jet if the cMVA discriminant is above a certain threshold, chosen such that the misidentification rate is about 1% for light-flavor quark and gluon jets, and about 13% for charm quark jets. The b jet tagging efficiency for this working point is about 66%.
The missing transverse momentum vector p miss T is computed as the negative vector sum of the transverse momenta of all the PF candidates in an event, and its magnitude is denoted as p miss T [54]. The p miss T is modified to account for corrections to the energy scale of the reconstructed jets in the event.
The candidate 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 [49,50] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets.
Muons are reconstructed as tracks in the muon system that are matched to the tracks reconstructed in the inner silicon tracking system [55]. The leading muon is required to have p T > 20 GeV, while the subleading muon must have p T > 15 (10) GeV in the bb νν (bb jj) channel. Muons are required to be reconstructed in the HLT fiducial volume, i.e., with |η| < 2.4, to ensure that the offline selection is at least as restrictive as the HLT requirements. The selected muons are required to satisfy a set of identification requirements based on the number of spatial measurements in the silicon tracker and in the muon system and the fit quality of the combined muon track [45], and are required to be consistent with originating from the primary vertex.
Electrons are reconstructed by matching tracks in the silicon tracker to the clusters of energy deposited in the ECAL [56]. The leading (subleading) electron is required to have p T > 25 (15) GeV and |η| < 2.5 to be within the geometrical acceptance, excluding candidates in the range 1.4442 < |η| < 1.5660, which is the transition region between the ECAL barrel and endcaps, because the reconstruction of an electron in this region is poor compared to other regions. Electrons are required to pass an identification requirement based on an MVA [57] technique that combines information from various observables related to the shower shape in the ECAL and the quality of the matching between the tracks and the associated ECAL clusters [56]. They are further required to be consistent with originating from the primary vertex. Candidates that are identified as originating from photon conversions in the material of the detector are removed.
Both muons and electrons have a requirement that the lepton relative isolation, defined in Eq.(1), be less than 0.25 (0.15) and 0.15 (0.06), respectively, for the bb jj (bb νν) channel. In Eq.(1), the sums run over charged hadrons originating from the primary vertex of the event, neutral hadrons, and photons inside a cone of radius ∆R = 3) around the direction of the muon (electron), where φ is the azimuthal angle in radians.
The isolation includes a correction for pileup effects, Corr PU . For electrons, Corr PU = ρA eff , where ρ is the average transverse momentum flow density, calculated using the jet area method [58], and A eff is the geometric area of the isolation cone times an η-dependent correction factor that accounts for residual pileup effects. For muons, Corr PU = 0.5 ∑ PU p T , where the sum runs over charged particles not associated with the primary vertex and the factor 0.5 corresponds to an approximate average ratio of neutral to charged particles in the isolation cone [59].
Simulated background and signal events are corrected with scale factors for differences observed between data and simulation, in trigger efficiencies, in lepton p T -and η-dependent identification and isolation efficiencies, and in b tagging efficiencies.

Event selection in the bb jj channel
After selection of the candidate physics objects, an initial event selection is performed by requiring at least two same-flavor leptons (muons or electrons) in each event. The two leptons are required to be oppositely charged. The invariant mass of the two leptons, m , is required to be larger than 15 GeV. Four of the jets in an event are designated as the h and Z boson decay products. These jets are required to have p T > 20 GeV and at least one of those must be b tagged with a minimum requirement on the b tagging discriminant, that is looser than the requirement in the final selection. We refer to this selection as the preselection.
Since the signal contains two b jets from the decay of a Higgs boson, and two jets of any flavor from the decay of a Z boson, it is important to carefully categorize the jets in the event. Starting from a collection of jets identified as described above, the information from the b tagging discriminant, as well as the kinematic properties of the jets, are taken into account when assigning jets as each particle's decay products.
The following selection is applied to identify the b jets originating from the decay of the Higgs boson. The two jets with the highest b tagging scores above a certain threshold are assigned to the decay of the Higgs boson. If only one jet is found that meets the minimum b tagging score value, a second jet that leads to an invariant mass closest to 125 GeV is selected. If no jets with b tagging scores above threshold are found, the two jets whose invariant mass is closest to 125 GeV are chosen.
After jets are assigned to the decay of h → bb, from the remaining jets the two jets with fourobject invariant mass M( jj) closest to 125 GeV are assigned to the decay of the Z boson.
After preselection, additional requirements are imposed. At least one of the four jets assigned as the decay products of the h or Z boson must satisfy the b tagging requirement, to increase the signal-to-background ratio. To impose orthogonality with the bb νν decay channel, upper limits on the p miss T are imposed as follows: p miss T < 40, 75, and 100 GeV for the m X of 260-350, 350-650, and ≥650 GeV, respectively. We refer to this selection as the final selection in the bb jj channel.
After the final selection, twenty-two variables that exploit the differences in kinematic and angular distributions between the signal and background processes are combined into a boosted decision tree (BDT) discriminant [60]. In the m X range of 260-300 GeV, the most important variables are m , the separation between the leading lepton and leading b tagged jet ∆R 1b1 , and the invariant mass of the pair of b tagged jets m h bb . In the m X range of 350-550 GeV, m h bb is the most important variable, while m becomes less important, and the separation between the pair of leptons ∆R gradually becomes more important when the m X increases. For the m X higher than 550 GeV, ∆R becomes the most important variable followed by m h bb and the separation between the pair of b tagged jets ∆R h bb . The BDTs are configured to use stochastic gradient boosting with the binomial log-likelihood loss function. The software package TMVA [57] is used for BDT implementation, training, and application.
The BDT is trained using all background processes described in Section 4, excluding the multijet background. In each lepton channel and for each spin hypothesis, one BDT is trained for each simulated signal m X . In the training, signal events include samples from the two neighboring mass points, in addition to the targeted mass point. In total, 48 BDTs are trained. These BDT distributions for data and expected backgrounds are used as the final discriminating variable in the analysis.

Background estimation in the bb jj channel
The main processes that can mimic the signature of the signal in the bb jj channel are Z/γ * +jets and tt, with smaller contributions from QCD multijets, diboson+jets, W+jets, and SM Higgs boson production.
The contribution from the principal background, Z/γ * +jets, is estimated with simulated events normalized to the data at the preselection level in the Z boson-enriched control region 80 < m < 100 GeV. The contribution from tt is estimated in a similar manner, with the tt-enriched control region defined by m > 100 GeV, and p miss T > 100 GeV. The data-to-simulation normalization factors derived from the two control regions are R Z = 1.14 ± 0.01 (stat) and R tt = 0.91 ± 0.01 (stat) in the muon channel and R Z = 1.24 ± 0.01 (stat) and R tt = 0.97 ± 0.02 (stat) in the electron channel. These normalization factors are found to be consistent between lepton flavors when applying lepton-specific systematic variations.
The contribution from QCD multijet processes is determined from data with a method that exploits the fact that neither signal events nor events from other backgrounds produce final states with same-sign leptons at any significant level. Data events with same-sign isolated leptons are used to model the shape of the multijet background, after all non-QCD sources of background contributing to this selection are subtracted using simulation. The yield in this region is normalized with the ratio of the number of events with nonisolated OS leptons to the number of events with nonisolated same-sign leptons. Here, nonisolated leptons are those muons (electrons) that fail the relative isolation requirements described in Section 5.1. All non-QCD sources of background, estimated with simulated events, are subtracted from the numerator and the denominator before computing the ratio.
The contributions from diboson+jets, W+jets, and SM Higgs boson production are estimated from simulation.

Event selection in the bb νν channel
Candidate events in the bb νν channel are reconstructed from the physics objects, as described above. The two leptons (muons or electrons) are required to have OS, and the invariant mass of the two leptons, m , is required to exceed 76 GeV. One of the Higgs bosons is formed from the pair of b jets with the highest output value of the b tagging discriminant, and the second Higgs boson is reconstructed as a combination of the two charged leptons and the p miss T , representing the visible and invisible decays products, respectively, of the pair of Z bosons. The requirement on m reduces the contribution from resonant X → hh production in the bbWW final state, and makes this measurement orthogonal to a previous bbWW search [18], where only events with m below 76 GeV were considered.
For the Higgs boson decaying to a pair of Z bosons, the two neutrinos are not reconstructed in the detector, and a pseudo invariant mass of the Higgs boson is used to approximate the incomplete momentum four-vector of the H. The pseudo invariant mass is formed from the momenta of the two charged leptons coming from one of the Z bosons and the four-vector (p miss T , p miss T ) approximating that of the two-neutrino system coming from the other of the Z bosons, where the z component of p miss T is zero. While the true invariant mass of the pair of neutrinos is not zero but is equal to the invariant mass of the parent Z boson, that boson is off the mass shell and has relatively low mass.
In order to suppress the backgrounds from the Z/γ * +jets and QCD multijet processes as well as from the SM Higgs boson production via the Zh process, a requirement is imposed on the minimum p miss T , which is 40 (75) GeV for the m X of 260-300 (350-600) GeV, and 100 GeV for higher m X .
Three regions, a signal region (SR) and two control regions (CRs), are further defined using m and the invariant mass m h bb of the two b jets. The SR is defined by the requirements 76 < m < 106 GeV and 90 < m h bb < 150 GeV. A first CR, dominated by tt events, is defined by m > 106 GeV and 90 < m h bb < 150 GeV. A second CR, enriched in Z/γ * +jets events, is defined by requiring 20 < m h bb < 90 GeV or m h bb > 150 GeV, and 76 < m < 106 GeV. The two CRs and the SR are used to estimate the backgrounds in the SR via a simultaneous fit.
To further differentiate signal from backgrounds in the SR, a BDT discriminant is trained using all simulated signal and background processes described in Section 4. Of the nine input distributions to the BDT, the most important variables in the low-mass range are the separation between the pair of b tagged jets ∆R h bb , p miss T , and m h bb . In the high-mass region, m h bb and ∆R h bb are also the most significant, together with the separation between the pair of charged leptons ∆R , which becomes more important as the resonance mass increases. Two BDTs are trained for each lepton channel and each resonance spin hypothesis, one for m X in the range of 250-450 GeV, and another one for the m X above 450 GeV. A minimum BDT value is required for candidates in the SR, optimized for each narrow m X hypothesis to yield the best 95% confidence level (CL) expected upper limit on resonant production. The BDTs are configured with the same classification and loss function parameters described in Section 5.2.
Finally, a quantity closely correlated with the energy-momentum four-vector of the hh system is constructed as the vector sum of the of the two leptons, two b jets, and the four-vector formed as (p miss T , p miss T ) for the neutrinos, as described above. Subsequently, the pseudo transverse mass of the hh system is defined as M T (hh) = √ E 2 − p 2 z , where E and p z are the energy and the z component of the combined four-vector.
The M T (hh) distributions for data and expected backgrounds, in the combined signal and CRs, will be used as the final discriminating variable in the analysis.
After the event selection in this channel is applied, the signal hh events in the SR come predominantly from the decays with the bbZZ intermediate state (80%) with a smaller contribution from the bbW + W − intermediate state (20%).

Background estimation in the bb νν channel
The dominant sources of background in the bb νν channel are tt and Z/γ * +jets production. Several other background processes contribute, including single top quark and diboson production, and SM Higgs boson production in association with a Z boson. While these are typically minor backgrounds, their contribution can vary over the m X range. The QCD multijet background is negligible across the full mass range because of the stringent selection on m .
The event yields in the signal and two CRs, which are dominated by tt and Z/γ * +jets events, are determined from data. The corresponding normalizations of the simulated M T (hh) distributions are free parameters in the simultaneous fit of all three regions. The remaining backgrounds are estimated from simulation and normalized according to their theoretical cross sections.

Systematic uncertainties
The dominant source of systematic uncertainty in this analysis is the jet energy scale (JES) uncertainty, which is of the order of a few percent and is estimated as a function of jet p T and η [52]. The η-dependent jet energy resolution (JER) correction factors are varied by ±1 standard deviation in order to estimate the effect of the uncertainty. Uncertainties in the JES are propagated to the calculation of p miss T . A residual p miss T uncertainty of 3% is applied in the bb νν channel to take into account the effect, at low p miss T , of the unclustered energy from neutral hadrons and photons that do not belong to any jet, and from jets with p T < 10 GeV.
An uncertainty of 2% per muon in the muon reconstruction, identification, and isolation requirements, as well as a 1% per muon uncertainty in the muon HLT efficiency are assigned [55]. A per-muon uncertainty due to measured differences of tracking efficiency in data and simulation is estimated to be 0.5% for muon p T < 300 GeV and 1.0% for muon p T > 300 GeV [61]. Per-electron uncertainties in the efficiency for electron trigger, identification and isolation requirements, estimated by varying the scale factors within their uncertainties, are applied. The uncertainties in the efficiency scale factors are generally <2% for trigger and <3% for identification and isolation [56]. The effect of the variations on the yield of the total background is <1%. Uncertainties in the data-to-simulation correction factors of the b tagging and of light flavor mis-tagging efficiencies are included.
Normalization and shape uncertainties are assigned to the modeling of the backgrounds. An uncertainty in the shape of the signal and background models is determined by varying the factorization and the renormalization scales between their nominal values and 0.5 to 2.0 times the nominal values in the simulated signal and background samples. The variations where one scale increases and the other decreases are not considered. Each of the remaining variations of the renormalization and the factorization scales are considered, and the maximum variation among all the samples with respect to the nominal sample used in the analysis is taken as the systematic uncertainty, which is found to be 5-7% depending on the process. An uncertainty in the signal acceptance and background acceptance and cross section due to PDF uncertainties and to the value chosen for the strong coupling constant is estimated by varying the NNPDF set of eigenvectors within their uncertainties, following the PDF4LHC prescription [62]. Statistical uncertainties in the simulated samples for Z/γ * +jets and tt background estimates result in uncertainties on the data-derived normalization factors in the bb jj channel.
An uncertainty of 2.5% is assigned to the determination of the integrated luminosity [63]. The uncertainty in the PU condition and modeling is assessed by varying the inelastic pp cross section from its central value by ±4.6% [64].
All the uncertainties discussed are applied to all background and signal simulated samples. The sensitivity of the presented search is limited by the statistical uncertainties.

Results
Results are obtained by performing a binned maximum likelihood fit of the BDT distributions for the bb jj channel, and of the hh pseudo transverse mass simultaneously in the SR and two CRs for the bb νν channel.
The data and background predictions at final selection level in the bb jj channel are shown in Fig. 1, for the distributions of the BDT discriminant for signal masses of 500 and 1000 GeV, in the muon and electron final states. Studies performed on all 48 BDT discriminants show stability of the trainings with no evidence of bias or overtraining.

CMS
Signal region channel ν ν bbee Figure 2: Pseudo transverse mass of the reconstructed hh candidates, in the bb νν channel, for data, simulated spin-2 RS1 graviton signal with a mass of 300 GeV, and simulated backgrounds scaled according to the fit results. The upper and lower rows correspond to the muon and electrons channels. For each row, the left and middle plots are for the Z/γ * +jets and tt control regions, and the right is for the signal region. The signals are normalized to 1 pb for the pp → X → hh process. The shaded area represents the combined statistical and systematic uncertainties in the background estimate. Figure 2 shows the hh pseudo transverse mass distributions in the data, background estimates, and spin-2 RS1 graviton for the 300 GeV mass hypothesis, after the final selection in the bb νν channel.
The systematic uncertainties are represented by nuisance parameters that are varied in the fit according to their probability density functions, prescribed as follows. A log-normal probability density function is assumed for the nuisance parameters affecting the event yields of the various background contributions, whereas systematic uncertainties that affect the distributions are represented by nuisance parameters whose variation is a vertical interpolation in each bin with a sixth-order polynomial for upward and downward shifts of one standard deviation, and linearly outside of that [65].
The statistical uncertainty from the limited number of events in the simulated samples is taken into account, for each bin of the discriminant distributions, by assigning a nuisance parameter to scale the sum of the process yields in that bin according to the statistical uncertainty using the Barlow-Beeston "lite" prescription [66,67].
In both channels the data distributions are well reproduced by the SM background processes. Upper limits on the resonance production cross section are set, using the asymptotic CL s modified frequentist approach [68][69][70].
The observed and expected 95% CL upper limits on σ(pp → X → HH → bbZZ) in the bb jj and bb νν channels as a function of m X are shown in Fig. 3, together with the NLO predictions for the RS1 radion, RS1 KK graviton, and N2HDM resonance production cross sections, where H can represent either the SM Higgs boson or an additional Higgs boson from an extended electroweak sector. As two different BDTs are defined for the search in the low-and high-mass ranges of the bb νν channel, the limit calculation is performed with both of the BDTs at the boundary of the two ranges, around 450 GeV, where a discontinuity is seen.
Combined 95% CL upper limits from both channels on σ(pp → X → HH → bbZZ) as a function of m X , are shown in Fig. 4, together with the theoretical predictions for the RS1 radion and RS1 KK graviton. In the m X range between 260 and 1000 GeV, limits on the production cross section times branching fraction of RS1 radion and RS1 KK graviton range from 0.1 to 5.0 and 0.1 to 3.6 pb, respectively. In the spin-0 case, the predictions of the N2HDM model with tan β = 0.5 and 2.0 are shown, for all H 3 → H 1 H 1 /H 1 H 2 /H 2 H 2 → bbZZ decays. In the tan β = 0.5 case, the model can be excluded with H 3 in the m X range of 360-620 GeV. In comparison to previous searches, we achieve a sensitivity to the RS1 radion and RS1 KK graviton models comparable to the difference in hh decay branching ratios.
Finally, the results are also interpreted as a function of both the m X and λ R (k) for the radion (graviton) case, with λ R < 0.3 TeV (k > 0.6) excluded for all of the m X considered, as shown in Fig. 5.

Summary
A search for the production of a narrow-width resonance decaying into a pair of Higgs bosons decaying into the bbZZ channel is presented. The analysis is based on data collected with the CMS detector during 2016, in proton-proton collisions at the LHC, corresponding to an integrated luminosity of 35.9 fb −1 . The final states considered are the ones where one of the Z bosons decays into a pair of muons or electrons, and the other Z boson decays either to a pair of quarks or a pair of neutrinos. Upper limits at 95% confidence level are placed on the production of narrow-width spin-0 or spin-2 particles decaying to a pair of Higgs bosons, in models with and without an extended Higgs sector. For a resonance mass range between 260 and 1000 GeV, limits on the production cross section times branching fraction of a spin-0 and spin-2 resonance range from 0.1 to 5.0 pb and 0.1 to 3.6 pb, respectively. These results set limits in parameter space in bulk Randall-Sundrum radion, Kaluza-Klein excitation of the graviton, and N2HDM models. For specific choices of parameters the N2HDM can be excluded in a mass range between 360 and 620 GeV for a resonance decaying to two Higgs bosons. This is the first search for Higgs boson resonant pair production in the bbZZ channel.     Figure 4: Expected (black dashed line) and observed (black solid line) limits on the cross section of resonant HH production times the branching fraction of HH → bbZZ as a function of the mass of the resonance for the combination of the bb jj and bb νν channels, where H can represent either the SM Higgs boson or an additional Higgs boson from an extended electroweak sector. The spin-0 case is shown on the left and the spin-2 case is shown on the right. The expected limits for each individual channel are shown with a red dashed line for the bb jj channel and blue dashed line for the bb νν channel. The red solid lines show the theoretical prediction for the cross section of an RS1 radion with λ R = 1 TeV and kL = 35 (left) and an RS1 KK graviton withk = 0.1 (right). In the spin-0 case only, the blue (green) line shows the decays of H 3 → H 1 H 1 /H 1 H 2 /H 2 H 2 → bbZZ in the N2HDM formulation, with tan β = 0.5 (2.0), the scalar H 3 vev set to 45 GeV, and the mixing angles α 1 , α 2 , α 3 set to 0.76, 0.48, and 1.00, respectively. The correction factor based on the relative partial width of H 3 to two gluons is around 3.0 (0.7) for tan β = 0.5 (2.0). The vertical black dashed line indicates the resonance mass of 450 GeV, a mass point where the BDT used in the analysis is switched from the one trained for low mass resonance to the one trained for high mass resonance.