Search for higgsino pair production in pp collisions at sqrt(s) = 13 TeV in final states with large missing transverse momentum and two Higgs bosons decaying via H to bb-bar

Results are reported from a search for new physics in 13 TeV proton-proton collisions in the final state with large missing transverse momentum and two Higgs bosons decaying via H to bb-bar. The search uses a data sample accumulated by the CMS experiment at the LHC in 2016, corresponding to an integrated luminosity of 35.9 inverse femtobarns. The search is motivated by models based on gauge-mediated supersymmetry breaking, which predict the electroweak production of a pair of higgsinos, each of which can decay via a cascade process to a Higgs boson and an undetected lightest supersymmetric particle. The observed event yields in the signal regions are consistent with the standard model background expectation obtained from control regions in data. Higgsinos in the mass range 230-770 GeV are excluded at 95% confidence level in the context of a simplified model for the production and decay of approximately degenerate higgsinos.

We perform a search for processes leading to Higgs boson pair production in association with large missing transverse momentum, p miss T . Each Higgs boson is reconstructed in its dominant decay mode, H → bb, which has a branching fraction of approximately 60%. Such a signature can arise, for example, in models based on SUSY, in which an electroweak process can lead to the production of two SUSY particles, each of which decays into a Higgs boson and another particle that interacts so weakly that it escapes detection in the apparatus. In this paper, we denote the particle in the search signature as H because it corresponds to the particle observed by ATLAS and CMS. However, in the context of SUSY models such as the minimal supersymmetric standard model (MSSM), this particle is usually assumed to correspond to the lighter of the two CP-even Higgs particles, denoted as h. The search uses an event sample of proton-proton (pp) collision data at √ s = 13 TeV, corresponding to an integrated luminosity of 35.9 fb −1 , collected by the CMS experiment at the CERN LHC. Searches for this and related decay scenarios have been performed by ATLAS [19,20] and CMS [15][16][17] using 7 and 8 TeV data. The analysis described here is based on an approach developed in Ref. [15].
While the observation of a Higgs boson completes the expected spectrum of SM particles, the low value of its mass raises fundamental questions that suggest the existence of new physics beyond the SM. Assuming that the Higgs boson is a fundamental (that is, noncomposite) spin-0 particle, stabilizing the electroweak scale (and the Higgs boson mass with it) is a major theoretical challenge, referred to as the gauge hierarchy problem [21][22][23][24][25][26]. Without invoking new physics, the Higgs boson mass would be pulled by quantum loop corrections to the cutoff scale of the theory, which can be as high as the Planck scale. Preventing such behavior requires an extreme degree of fine tuning of the theoretical parameters. Alternatively, stabilization of the Higgs boson mass can be achieved through a variety of mechanisms that introduce new physics at the TeV scale, such as SUSY.
The class of so-called natural SUSY models [27][28][29][30][31] contains the ingredients necessary to stabilize the electroweak scale. These models are the object of an intensive program of searches at the LHC. In any SUSY model, additional particles are introduced, such that all fermionic (bosonic) degrees of freedom in the SM are paired with corresponding bosonic (fermionic) degrees of freedom in the extended theory. In natural SUSY models, certain classes of partner particles are expected to be light. These include the four higgsinos ( H 0 1,2 , H ± ); both top squarks, t L and t R , which have the same electroweak couplings as the left-(L) and right-(R) handed top quarks, respectively; the bottom squark with left-handed couplings ( b L ); and the gluino ( g). Of these, the higgsinos are generically expected to be the lightest. Furthermore, in natural scenarios, the four higgsinos are approximately degenerate in mass, so that transitions among these SUSY partners would typically produce only very soft (i.e., low momentum) additional particles, which would not contribute to a distinctive experimental signature.
In general, the gaugino and higgsino fields can mix, leading to mass eigenstates that are classified either as neutralinos ( χ 0 i , i = 1-4) or charginos ( χ ± i , i = 1-2). If the χ 0 1 is the lightest SUSY particle (LSP), it is stable in models that conserve R-parity [32] and, because of its weak  Figure 1: Diagram for the gauge-mediated symmetry breaking signal model, χ 0 1 χ 0 1 → HH G G (TChiHH), where G is a goldstino. The NLSPs χ 0 1 are not directly pair produced, but are instead produced in the cascade decays of several different combinations of neutralinos and charginos, as described in the text.
interactions, would escape experimental detection. Achieving sensitivity to scenarios in which the higgsino sector is nearly mass degenerate and contains the LSP poses a significant experimental challenge because the events are characterized by low-p T SM decay products and small values of p miss T [33][34][35]. Searches based on signatures involving initial-state radiation (ISR) recoiling against the pair produced higgsinos have already excluded limited regions of phase space for these scenarios [36][37][38][39][40]. However, achieving broad sensitivity based on this strategy is expected to require the large data samples that will be accumulated by the HL-LHC [41].
An alternative scenario arises, however, if the lightest higgsino/neutralino is not the LSP, but the next-to-lightest SUSY particle (NLSP) [42]. The LSP can be another particle that is generic in SUSY models, the goldstino ( G). The goldstino is the Nambu-Goldstone particle associated with the spontaneous breaking of global supersymmetry and is a fermion. In a broad range of models in which SUSY breaking is mediated at a low scale, such as gauge-mediated supersymmetry breaking (GMSB) models [43,44], the goldstino is nearly massless on the scale of the other particles and becomes the LSP. If SUSY is promoted to a local symmetry, as is required for the full theory to include gravity, the goldstino is "eaten" by the SUSY partner of the graviton, the gravitino (J = 3/2), and provides two of its four degrees of freedom. In the region of parameter space involving prompt decays to the gravitino, only the degrees of freedom associated with the goldstino have sufficiently large couplings to be relevant, so it is common to denote the LSP in either case as a goldstino. In these GMSB models, the goldstino mass is generically at the eV scale.
If the lighter neutralinos and charginos are dominated by their higgsino content and are thus nearly mass degenerate, their cascade decays can all lead to the production of the lightest neutralino, χ 0 1 (now taken to be the NLSP), and soft particles. Integrating over the contributions from various allowed combinations of produced charginos and neutralinos ( χ 0 therefore leads to an effective rate for χ 0 1 χ 0 1 production [45,46] that is significantly larger than that for any of the individual primary pairs, resulting in a boost to the experimental sensitivity. The higgsino-like NLSP would then decay via χ 0 1 → γ G, χ 0 1 → H G, or χ 0 1 → Z G, where the goldstino can lead to large p miss T . The branching fractions for these decay modes vary depending on a number of parameters including tan β, the ratio of the Higgs vacuum expectation values, and the branching fraction for χ 0 1 → H G can be substantial. As a consequence, the signature HH+p miss T with H → bb can provide sensitivity to the existence of a higgsino sector in the important class of scenarios in which the LSP mass lies below the higgsino masses. This article presents a search for higgsinos in events with p miss T > 150 GeV and at least three jets identified as originating from b quark hadronization (b-tagged jets). In each event, we reconstruct two Higgs boson candidates and define search regions within a Higgs boson mass window. The background is dominated by tt production at low and intermediate p miss T , and by Z → νν production in association with b quarks at high p miss T . The background is estimated entirely from data control regions corresponding to events with two b-tagged jets and events with three or four b-tagged jets outside the Higgs boson mass window. Systematic uncertainties on the background prediction are derived from both dedicated data control regions for tt, Z → νν and QCD multijet production as well as from the simulation of the background events in the search regions.
Results are interpreted in the simplified model framework [47][48][49] using the model shown in Fig. 1, hereafter referred to as TChiHH. In this model, two χ 0 1 NLSPs are produced, each decaying via χ 0 1 → H G. The cross-section is calculated under the assumption that at least one of the χ 0 1 NLSPs is produced in a cascade decay of χ 0 2 or χ ± 1 , as described above. This situation arises when the mass splittings among charginos and neutralinos are large enough that the decays to χ 0 1 occur promptly [50], while also small enough that the additional soft particles fall out of acceptance.

The CMS detector
The central feature of the CMS detector is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are the tracking and calorimeter systems. The tracking system, composed of silicon-pixel and silicon-strip detectors, measures charged particle trajectories within the pseudorapidity range |η| < 2.5, where η ≡ − ln[tan(θ/2)] and θ is the polar angle of the trajectory of the particle with respect to the counterclockwise proton beam direction. A lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections, provide energy measurements up to |η| = 3. Forward calorimeters extend the pseudorapidity coverage provided by the barrel and endcap detectors up to |η| = 5. Muons are identified and measured within the range |η| < 2.4 by gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. The detector is nearly hermetic, permitting the accurate measurement of p miss T . A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, is given in Ref. [51].

Simulated event samples
Several simulated event samples are used for modeling the SM background and signal processes. While the background estimation in the analysis is performed from control samples in the data, simulated event samples are used to estimate uncertainties, as well as to build an understanding of the characteristics of the selected background events.
The production of tt+jets, W+jets, Z+jets, and quantum chromodynamics (QCD) multijet events is simulated with the Monte Carlo (MC) generator MADGRAPH5 aMC@NLO 2.2.2 [52] in leadingorder (LO) mode [53]. Single top quark events are modeled with POWHEG v2 [54,55] for the tchannel and tW production, and MADGRAPH5 aMC@NLO at next-to-leading order (NLO) [56] for the s-channel. Additional small backgrounds, such as tt production in association with bosons, diboson processes, and tttt, are also produced at NLO with either MADGRAPH5 aMC@NLO or POWHEG. All events are generated using the NNPDF 3.0 [57] set of parton distribution functions. Parton showering and fragmentation are performed with the PYTHIA 8.205 [58] generator with the underlying event model based on the CUETP8M1 tune [59]. The detector simulation is performed with GEANT4 [60][61][62]. The cross sections used to scale simulated event yields are based on the highest order calculation available [54,55,[63][64][65][66][67][68][69][70][71], which for the most part correspond to NLO or next-to-NLO precision.
Signal events for the TChiHH simplified model are generated for 36 values of the higgsino mass between 127 and 1000 GeV. The mass points are denoted as TChiHH(m χ 0 1 ,m G ), where m χ 0 1 is the higgsino mass and m G is the mass of the LSP, both in units of GeV. While the value of m G is fixed to 1 GeV in the simulation for technical reasons, the resulting event kinematics are consistent with an approximately massless LSP such as the goldstino in GMSB. The yields are normalized to the NLO + next-to-leading logarithmic (NLL) cross section [45,46]. The production cross sections are calculated in the limit of mass degeneracy among higgsinos, χ 0 1 , χ 0 2 , and χ ± 1 . All other supersymmetric partners of the SM particles are assumed to be heavy (100 TeV) and thus essentially decoupled, a simplification that has a very small impact on higgsino pair production (e.g., the cross section changes less than 2% when the masses of the gluino and squarks are lowered down to 500 GeV). Following the convention of real mixing matrices and signed neutralino masses [72], we set the sign of the mass of χ 0 1 ( χ 0 2 ) to +1 (−1). The lightest two neutralino states are defined as symmetric (anti-symmetric) combinations of higgsino states by setting the product of the elements N i3 and N i4 of the neutralino mixing matrix N to +0.5 (−0.5) for i = 1 (2). The elements U 12 and V 12 of the chargino mixing matrices U and V are set to 1. All chargino and neutralino decays in the simplified model are taken to be prompt, although the lifetimes of particles in a physical model would depend on the mass splitting between the higgsino-like states, which become long-lived in the limit of degenerate masses. Both Higgs bosons in each event are forced to decay to bb, which is accounted for by scaling the signal event yields with the branching fraction of 58.24% [73]. The signal contribution from Higgs boson decays other than H → bb is small in this analysis and is ignored. Signal events are generated in a manner similar to that for the SM backgrounds, with the MAD-GRAPH5 aMC@NLO 2.2.2 generator in LO mode using the NNPDF 3.0 set of parton distribution functions and followed by PYTHIA 8.205 for showering and fragmentation. The detector simulation is performed with the CMS fast simulation package [74] with scale factors applied to compensate for any differences with respect to the full simulation. Finally, to model the presence of additional pp collisions from the same beam crossing as the primary hard-scattering process or another beam crossing adjacent to it ("pileup" interactions), the simulated events are overlaid with multiple minimum-bias events, which are also generated with the PYTHIA 8.205 generator with the underlying event model based on the CUETP8M1 tune.

Object definitions
The reconstruction of physics objects in an event proceeds from the candidate particles identified by the particle-flow (PF) algorithm [75], which uses information from the tracker, calorimeters, and muon systems to identify the candidates as charged or neutral hadrons, photons, electrons, or muons. The reconstructed vertex with the largest value of summed physics-object p 2 T , with p T denoting transverse momentum, is taken to be the primary pp interaction vertex (PV). The physics objects used in this context are the objects returned by a jet finding algorithm [76,77] applied to all charged tracks associated with the vertex under consideration, plus the corresponding associated p miss T .
The charged PF candidates associated with the PV and the neutral PF candidates are clustered into jets using the anti-k T algorithm [76] with a distance parameter R = 0.4, as implemented in the FASTJET package [77]. The jet momentum is determined as the vectorial sum of all particle momenta in the jet. Jet energy corrections are derived based on a combination of simulation studies, accounting for the nonlinear detector response and the presence of pileup, together with in situ measurements of the p T balance in dijet and γ+jet events [78]. The resulting calibrated jet is required to satisfy p T > 30 GeV and |η| ≤ 2.4. Additional selection criteria are applied to each event to remove spurious jet-like features originating from isolated noise in certain HCAL regions [79].
A subset of the jets are "tagged" as originating from b quarks using DEEPCSV [80], a new b tagging algorithm based on a deep neural network with four hidden layers [81]. We use all three of the DEEPCSV algorithm working points: loose, medium, and tight, defined by the values of the discriminator requirement for which the rates for misidentifying a light-flavor jet as a b jet are 10%, 1%, and 0.1%, respectively. The b tagging efficiency for jets with p T in the 80-150 GeV range is approximately 86%, 69%, and 51% for the loose, medium, and tight working points, respectively, and gradually decreases for lower and higher jet p T . The simulation is reweighted to compensate for any differences with respect to data based on measurements of the b tagging efficiency and mistag rate for each working point in dedicated data samples [80].
The missing transverse momentum, p miss T , is given by the magnitude of p miss T , the negative vector sum of the transverse momenta of all PF candidates [82,83], adjusted for known detector effects by taking into account the jet energy corrections. Filters are applied to reject events with well defined anomalous sources of p miss T arising from calorimeter noise, dead calorimeter cells, beam halo, and other effects.
Since the targeted signature is fully hadronic, contamination from final states involving leptons in the search region is suppressed by vetoing events with reconstructed lepton candidates. Electrons are identified by associating a charged particle track with an ECAL supercluster [84] and are required to have p T > 10 GeV and |η| < 2.5. Muons are identified by associating tracks in the muon system with those found in the silicon tracker [85] and are required to satisfy p T > 10 GeV and |η| < 2.4.
To preferentially consider only leptons that originate in the decay of W and Z bosons, leptons are required to be isolated from other PF candidates using an optimized version of the "miniisolation" [86,87]. The isolation I mini is obtained by summing the transverse momentum of the charged hadrons, neutral hadrons, and photons within ∆R ≡ √ (∆φ) 2 + (∆η) 2 < R 0 of the lepton momentum vector p , where φ is the azimuthal angle in radians and R 0 is given by 0.2 for p T ≤ 50 GeV, (10 GeV)/p T for 50 < p T < 200 GeV, and 0.05 for p T ≥ 200 GeV. Electrons (muons) are then required to satisfy I mini /p T < 0.2 (0.1).
As described in Section 5, the dominant background arises from the production of single-lepton tt events in which the lepton is a τ decaying hadronically, or is a light lepton that is either not reconstructed or fails the lepton selection criteria, including the p T threshold and the isolation requirements. To reduce this background, we veto events with any additional isolated tracks corresponding to leptonic or hadronic PF candidates. To reduce the influence of tracks originating from pileup, isolated tracks are considered only if their closest distance of approach along the beam axis to a reconstructed vertex is smaller for the primary event vertex than for any other vertex.
The requirements for the definition of an isolated track differ slightly depending on whether the track is identified as leptonic or hadronic by the PF algorithm. For leptonic tracks, we require p T > 5 GeV and I trk < 0.2, where I trk is the scalar p T sum of other charged tracks within ∆R < 0.3 of the primary track, divided by the p T value of the primary track. For hadronic tracks, we apply slightly tighter requirements to reduce hadronic (non-τ) signal loss: p T > 10 GeV and I trk < 0.1. To minimize the signal inefficiency due to this veto, isolated tracks are considered only if they are consistent with originating from a W boson decay, specifically, if the transverse mass of the track and the missing transverse momentum satisfy where p trk T is the transverse momentum of the track and ∆φ p trk T , p miss T is the azimuthal separation between the track and p miss T .
The majority of QCD multijet events containing high p miss T have at least one jet with undermeasured momentum and thus a spurious momentum imbalance. A signature of such an event is a jet closely aligned in direction with the p miss T vector. To suppress this background, we place the following requirements on the angle ∆φ i between the ith highest p T jet and p miss T for i = 1, 2, 3, 4: ∆φ 1 > 0.5, ∆φ 2 > 0.5, ∆φ 3 > 0.3, and ∆φ 4 > 0.3. These conditions are hereafter collectively referred to as the high ∆φ requirement.
The number of jets satisfying the criteria described above is denoted as N jets , while the numbers of these jets tagged with the loose, medium, and tight b tagging working points are labeled as N b,L , N b,M , and N b,T , respectively. By definition, the jets identified by each b tagging working point form a subset of those satisfying the requirements of looser working points.

Definition of b tag categories
To optimize signal efficiency and background rejection, we define the following mutually exclusive b tag categories: The 2b category is used as a control sample to determine the kinematic shape of the background. Most of the signal events lie in the 3b and 4b categories. This categorization is found to have superior performance with respect to other combinations of working points. For instance, the simpler option of only using medium b tags results in a 2b control sample that has a larger contribution from QCD multijet production and a 4b sample with smaller signal efficiency.
To study various sources of background with higher statistical precision, we also define the following b tag categories with looser requirements: We will use N b as a shorthand when discussing b tag categories as an analysis variable, and N b,L , N b,M , and N b,T when discussing numbers of b-tagged jets for specific working points.

Higgs boson pair reconstruction
The principal visible decay products in signal events are the four b jets from the decay of the two Higgs bosons. Additional jets may arise from initial-or final-state radiation as well as pileup. To reconstruct both Higgs bosons, we choose the four jets with the largest DEEPCSV discriminator values, i.e., the four most b-quark-like jets. These four jets can be grouped into three different pairs of Higgs boson candidates. Of the three possibilities, we choose the one with the smallest mass difference ∆m between the two Higgs boson candidate masses m H 1 and This method exploits the fact that signal events contain two particles of identical mass, without using the known value of the Higgs boson mass itself. Methods that use the known mass to select the best candidate tend to sculpt an artificial peak in the background.
Only events where the masses of the two Higgs boson candidates are similar, ∆m < 40 GeV, are kept. We then calculate the average mass as As discussed in Section 6, the search is then performed within the Higgs boson mass window defined as 100 < m ≤ 140 GeV.
After selecting the two Higgs boson candidates, we compute the distance ∆R between the two jets in each of the H → bb candidate decays. We then define ∆R max as the larger of these two values, In the typical configuration of signal events satisfying the baseline requirements, ∆R max is small because the Higgs bosons tend to have nonzero transverse boost and, thus, the two jets from the decay of a Higgs boson tend to lie near each other in η and φ. In contrast, for semileptonic tt background events, three of the jets typically arise from a top quark that decays via a hadronically decaying W boson while the fourth jet arises from a b quark from the other top quark decay. Therefore, three of the jets tend to lie within the same hemisphere while the fourth jet is in the opposite hemisphere. One of the Higgs boson candidates is thus formed from jets in both hemispheres, and ∆R max tends to be larger than it is for signal events.

Trigger and event selection
The data sample was obtained with triggers that require the online p miss T value to be greater than 100 to 120 GeV, the applied threshold varying with the instantaneous luminosity delivered by the LHC. This variable is computed with trigger-level quantities, and therefore has somewhat worse resolution than the corresponding offline variable. The trigger efficiency measured as a function of offline p miss T , in samples triggered by a high-p T isolated electron, rises rapidly from about 60% at p miss T = 150 GeV to 92% for p miss T = 200 GeV and to over 99% for p miss T > 300 GeV. The systematic uncertainty in the trigger efficiency is obtained by comparing the nominal efficiency with that found in different kinematic regions, with various reference triggers, and with the simulation. This uncertainty is about 7% for p miss T = 150 GeV and decreases to 0.7% for p miss T > 300 GeV.
Several data control samples are employed to validate the analysis techniques and to estimate systematic uncertainties in the background estimates. The control sample for the principal background from tt events requires exactly one electron or one muon, while the Z → νν background is studied with a control sample requiring two leptons consistent with a Z → + − decay. These data samples were obtained with triggers that require at least one electron or muon with p T greater than 27 or 24 GeV, respectively.

Trigger and event selection
Signal events have four b jets from the decay of two Higgs bosons and no isolated leptons, with any additional hadronic activity coming from initial-or final-state radiation. Thus, we select events with four or five jets, no leptons or isolated tracks, N b,T ≥ 2, p miss T > 150 GeV, high ∆φ, ∆m < 40 GeV, and ∆R max < 2.2. These selection requirements, listed in the top half of Table 1, are referred to as the baseline selection, while the bottom half of that table shows the further reduction in background in increasingly more sensitive search bins. The distributions of ∆m, ∆R max , and m in the 4b category are shown in Fig. 2 in data and simulation. The actual background prediction, however, is based on data control samples, as described in the next section. Table 1: Event yields obtained from simulated event samples scaled to an integrated luminosity of 35.9 fb −1 , as the event selection criteria are applied. The category "tt+X" is dominated by tt (98.5%), but also includes small contributions from tttt, ttW, ttZ, ttH, and ttγ backgrounds. The category "V+jets" includes Z+jets and W+jets backgrounds in all their decay modes. The category "Other" includes ZZ, WZ, WW, WH(→ bb), and ZH(→ bb) processes. The event selection requirements listed up to and including ∆R max < 2.2 are defined as the baseline selection. The trigger efficiency is applied as an event weight and is first taken into account in the p miss T > 150 GeV row. The uncertainties in the "Total SM bkg." column is statistical only. The columns corresponding to the yields for three signal benchmark points are labeled by TChiHH(m χ 0 1 ,m G ), with m χ 0 1 and m G in units of GeV. The simulated samples for TChiHH(225,1), TChiHH(400,1), and TChiHH(700,1) are equivalent to 10, 100, and over 1000 times the data sample, respectively, so the statistical uncertainties in the signal yields are small.
Based on the simulation, after the baseline selection, more than 85% of the remaining SM background arises from semileptonic tt production. Approximately half of this contribution corresponds to tt events with an electron or a muon in the final state that is either out of acceptance or not identified, while the other half involves final states with a hadronically decaying τ lepton. The contributions from events with a W or Z boson in association with jets (V+jets) are about 10% and are dominated by Z → νν decays. The background from QCD multijet events after the baseline selection is very small due to the combination of p miss T , ∆φ, and N b requirements. As shown in Fig. 2, the p miss T distribution of the signal is highly dependent on the higgsino mass. To further enhance the sensitivity of the analysis, we therefore subdivide the search region into four p miss   6 Background estimation 6

.1 Method
The background estimation method is based on the observation that the m distribution is approximately uncorrelated with the number of b tags. As shown in Fig. 3, the m shapes for the three b tag categories agree within the statistical uncertainty in the simulated samples. This behavior can be understood by noting that the background in all three b tag categories is dominated by events containing only two b quarks, with the additional b-tagged jets in the 3b and 4b categories being mistagged light-flavor or gluon jets. The background simulation indicates that only 20% (37%) of the events in the 3b ( 4b) selection have more than two b quarks. As a result, the four jets used to construct m in the 3b and 4b categories arise largely from the same fundamental processes as those with two b-tagged jets, and thus the shape of the average mass distribution is independent of N b for N b ≥ 2.
Taking advantage of this observation, we estimate the background contribution to each signal bin with an ABCD method [87] that employs m and the b tag categories as the two ABCD variables similarly to the 8 TeV analysis [15]. We define the Higgs boson mass window (HIG region) as the events with m within 100 to 140 GeV, and the Higgs boson mass sideband (SBD region) as the events with 0 < m < 100 GeV or 140 < m < 200 GeV. The mass window is chosen to optimize the signal sensitivity, taking into account the background distribution and the asymmetry in the Higgs boson mass resolution. The 3b and 4b SBD regions, together with the shape of the m distribution in the 2b category, are then used to determine the background in the signal-enriched 3b and 4b HIG regions independently for each p miss Here, µ bkg SBD,nb and µ bkg HIG,nb are the background rates for each b tag category (n = 2, 3, 4) in the SBD and HIG search regions, respectively, and R is the ratio of the background rate in the HIG κ 0 0  region to that in the SBD region. In the limit that the b tag category and m are uncorrelated, R is the same for the three b tag categories: The closure of the background estimation method, that is, the ability of Eq. (5) These κ factors measure the impact of any residual correlation between the b tag category and m . Figure 4 shows that the κ factors in simulation are consistent with unity for both the 3b and 4b regions across the full p miss T range, demonstrating the validity of the fundamental assumption of the ABCD method. In Section 7, we study the closure of the method in data control samples and estimate the associated systematic uncertainties in the background prediction.

Implementation
The method outlined in Section 6.1 is implemented with a likelihood function that incorporates the statistical and systematic uncertainties associated with the background prediction and the signal model, and also accounts for signal contamination in all control regions.
The terms in the likelihood function corresponding to the observed yields in all analysis regions, reflecting the parameterization of the ABCD method and the signal contributions to each bin, can be written as the following product of Poisson probability density functions (pdfs): Here, the index m runs over the four p miss T bins, the index n runs over the three b tag categories, N data are the observed data yields, µ sig are the expected signal rates, and r is the parameter quantifying the signal strength relative to the expected yield across all analysis bins. The four main floating parameters describing the fitted background for each p miss T bin m are the three sideband background rates µ bkg SBD,nb,m and the ratio R. The full likelihood function is given by the product of L ABCD , Poisson pdfs constraining the signal shape and its statistical uncertainty in each bin, and log-normal pdfs constraining nuisance parameters that account for the systematic uncertainties in the closure and the signal efficiency. These nuisance parameters were omitted from Eq. (8) for simplicity.
Following the approach in Ref. [87], the likelihood function is employed in two types of fits: the predictive fit, which allows us to more easily establish the agreement of the background predictions and the observations in the background-only hypothesis, and the global fit, which enables us to estimate the signal strength.
The predictive fit is realized by removing the terms of the likelihood corresponding to the observed yields in the signal regions, (HIG, 3b) and (HIG, 4b), and fixing the signal strength r to 0. As a result, we obtain a system of equations with an equal number of unknowns and constraints. For each p miss T bin, the four main floating parameters µ bkg SBD,2b , µ bkg SBD,3b , µ bkg SBD,4b , and R are determined by the four observations N data SBD,2b , N data SBD,3b , N data SBD,4b , and N data HIG,2b . Since the extra floating parameters corresponding to the systematic uncertainties are constrained by their respective log-normal pdfs, they do not contribute as additional degrees of freedom. The predictive fit thus converges to the standard ABCD method, and the likelihood maximization machinery becomes just a convenient way to solve the system of equations and to propagate the various uncertainties.
Conversely, the global fit includes the observations in the signal regions. Since in this case there are six observations and four floating background parameters in each p miss T bin, there are enough constraints for the signal strength r to be determined in the fit. The global fit also properly accounts for the signal yields in the control regions, using the signal shape across control and signal regions from the simulation.

Systematic uncertainties in the background prediction
The background estimation procedure described in Section 6 relies on the approximate independence of the m and N b distributions. In Sections 7.1, 7.2, and 7.3 we study this assumption for individual background processes in data and simulation by defining dedicated control regions for tt, Z+jets, and QCD multijet production. The overall level of closure in these control samples, better than 13% in all cases, is assigned as a systematic uncertainty for each of the main background sources separately. Additionally, these samples validate the closure in the simulation as a function of p miss T . If the background estimation method is valid for each separate background contribution, then it would also be valid for the full background composition as long as the relative abundance of each background component is independent of N b . In Section 7.4, we use these data control samples to quantify the validity of the simulation prediction that the background composition is independent of N b in each p miss T bin by examining the modeling of the p miss T and N b distributions for each background source.
Finally, in Section 7.5, we describe the prescription for assigning the total systematic uncer-tainty in the background prediction binned in p miss T and N b , taking into account both the findings from the data control sample studies and the closure of the method in the simulation. The latter is the dominant systematic uncertainty in this analysis.

Single-lepton tt control sample
To test whether the background estimation method works for tt events, we define a singlelepton control sample, which, like the search region, is dominated by single-lepton tt events and represents a similar kinematic phase space. Because the lepton is a spectator object as far as the ABCD method is concerned-it is neither involved in the construction of the m variable, nor correlated with the presence of additional b tags-this control sample should accurately capture any potential mismodeling of the m -N b correlation that may be present in the signal region. While this control region does not specifically probe semileptonic tt events involving a hadronically decaying τ lepton, the simulation shows that their m distribution in the signal region is the same as that of semileptonic tt events involving light leptons. This is expected because in most cases the τ lepton in these events is either out of acceptance or reconstructed as the jet with the smallest b-discriminator value and, as a result, it does not enter the m calculation.
For each of the four p miss T search bins, we construct a corresponding ABCD test in a singlelepton control sample, defined by the same selection requirements except for removing the lepton and the isolated track vetoes and instead requiring exactly one lepton with p T > 30 GeV (to reach trigger efficiency plateau) and m T ( p T , p miss T ) < 100 GeV (to avoid poorly reconstructed events). Given that the contamination from QCD multijet production in the single-lepton region is small, the ∆φ requirement is also removed to further improve the statistical power of the control sample. Since the lepton provides a way to trigger on events with lower p miss T , we add two additional p miss T bins, p miss T < 75 GeV and 75 < p miss T ≤ 150 GeV, allowing us to study the p miss T dependence of the closure in a wider range. In this control region, tt production accounts for over 90% of the events, except for the two highest p miss T bins, where the total contribution of single top quark and V+jets production can be as high as ∼25%. Figure 5 (left) shows the comparison of the m shapes between the data and the simulation.
As described in Section 6, since the 3b and 4b categories are dominated by events with two true B hadrons and one or two additional mistagged jets, similar jet topologies contribute to all b tag categories and thus the m distributions of the reconstructed b tag categories converge. We validate this assertion in the single-lepton control sample by examining the value of the κ factors. Figure 6 shows the overall closure of the method across bins of p miss T , both in the simulation and in data. We observe agreement within the statistical uncertainties, with κ values being consistent with unity across the full p miss T range for both data and simulation. This observation is also confirmed with larger statistical precision by repeating the test in a more inclusive sample obtained by removing the ∆R max requirement.
An overall uncertainty in the validity of the method in tt-like events is assigned based on the larger of the nonclosure and the statistical uncertainty in the closure test in data performed after integrating over the full p miss T range. The results, shown to the right of the solid line in Fig. 6, correspond to an uncertainty of 3% and 6% in the 3b and 4b bins, respectively.

Dilepton Z+jets control sample
As shown in Table 1, the second-largest background is Z+jets, with the Z boson decaying via Z → νν. Similarly to the tt case, we can validate the background estimation method for Z+jets events by constructing a closure test in a representative data control sample rich in Z → + −    decays. However, given the small branching fraction of Z → + − decays and the large tt contamination associated with a high-N b selection, we validate the method by constructing ABCD tests at lower b tag requirements, namely 1b/0b and 2b/1b. The Z → + − control sample is constructed in a similar manner to the search region. Events with 4 or 5 jets are selected, and the reconstruction of a pair of Higgs bosons proceeds as described in Section 4. We require two opposite-charge same-flavor signal leptons in the Z boson mass window, 80 < m( + − ) ≤ 100 GeV, with the p T of the leading and subleading lepton required to be greater than 40 and 20 GeV, respectively. We remove the lepton and isolated track vetoes and, since the dilepton selection makes the contamination from QCD multijet events negligible, we remove the ∆φ requirement. Since we do not expect genuine p miss T in Z → + − events, we additionally require p miss T < 50 GeV, which reduces the contamination from other processes from 20% to 10%.
We divide the sample in bins of p T ( + − ), ensuring kinematic correspondence with the Z → νν decays present in the various p miss T bins employed in the search region. Similarly to the singlelepton sample, the presence of leptons allows us to extend the closure test to lower values of p T ( + − ). Figure 5 (right) shows both the high purity of the sample and the excellent data-tosimulation agreement in the m shape.
The validity of the extrapolation of the method to a sample consisting of lower b tag multiplicities is supported by the observation that all jets in Z+jets events come from ISR, and thus their kinematic properties are largely independent of the flavor content of the event. This expectation is confirmed in data by examining the overall closure of the method in bins of p T ( + − ) for the 1b/0b and 2b/1b ABCD tests. The 1b/0b test, which has greater statistical power compared to the 2b/1b test and thus allows a better examination of any potential trends as a function of p T ( + − ), is shown in Fig. 7 for illustration.
Since we do not observe that the closure of the method has any dependence on p T ( + − ), we proceed to combine all the p T ( + − ) bins into one bin and repeat the closure test with improved statistical precision. In the 1b/0b ABCD test, we observe a statistically significant nonclosure of 11%, which may be due to higher order effects beyond the precision of this search. The 2b/1b ABCD test shows good closure but with a higher statistical uncertainty of 19%. We assign the larger uncertainty of 19% as the systematic uncertainty in the closure of the background estimate method for Z+jets events. The robustness of this result is further corroborated by similar checks in a more inclusive selection without the ∆R max requirement.

Low ∆φ QCD multijet control sample
Finally, to examine the validity of the ABCD method for QCD multijet events, we define a control region enriched with such events by inverting the ∆φ requirement. The high b tag multiplicity region of this control sample has a limited event yield and high tt contamination. To overcome these limitations, we exploit the fact that QCD multijet events, like Z+jets events, have similar kinematic properties regardless of their flavor content. We thus check the m -N b independence in lower b tag multiplicity regions by constructing the 1b/0b and 2b/1b ABCD tests. We observe good agreement between the data and the simulation for all p miss T bins. The maximum measured deviation of κ from unity in the inclusive bins equals 13%, which we assign as the systematic uncertainty in the closure of the background estimation method for QCD multijet events.

Impact of the background composition
Having evaluated the closure of the method for each individual background, we proceed to study the impact of mismodeling the relative abundance of the different background sources.
Since the m shape varies among background types, as shown for tt and Z+jets in Fig. 5, significant differences in the process admixture in the 2b category with respect to the 3b or 4b category will result in m -N b correlation and lead to the nonclosure of the method. From simulation, the background composition is expected to be independent of the b tag category. The validity of this prediction relies on the ability of the simulation to model the shape of the b tag category and p miss T distributions equally well for each background contribution.
From comparisons in the respective control samples, we indeed observe that the N b distribution for each of tt, Z+jets, and QCD multijet production is similarly well modeled by the simulation. The p miss T distribution in simulation is found to overestimate the data for large values of p miss T for tt and Z+jets events, while the opposite is observed for QCD multijet events. To provide an estimate of the impact of mismodeling the background composition, we reweight the simulation based on the data-to-simulation comparisons and then calculate the κ factors with the reweighted simulation, assigning 100% of the shift with respect to the nominal values as the uncertainty in the modeling of the background composition. The resulting uncertainty is found to be at most 4%.

Total systematic uncertainty determination
Based on the data control sample studies described in Sections 7.1-7.4, we assign a set of systematic uncertainties in the background prediction for each search bin as follows: 1. The closure uncertainty for each background process obtained in data control regions is propagated to the background predictions by varying the closure of the particular background in simulation in bins of p miss T and N b . The resulting shifts on the predictions, ranging from 1% to 10% increasing with p miss T , are assigned as systematic uncertainties with a 100% bin-to-bin correlation.
2. The level of nonclosure due to the relative abundance of each background component as a function of N b is estimated by comparing the change in κ in simulation before and after correcting the N b and p miss T distributions of each background source according to measurements in the data control samples. Its overall impact is 1-4% and it is taken as 100% correlated across the different analysis bins.
3. The closure studies in the data control samples with more inclusive selections show no evidence of p miss T dependence, but have insufficient statistical power at high p miss T using the default selection. Given this limitation and the extensive validation of the simulation in all control samples, we assign the larger of the statistical uncertainty and the nonclosure for each bin in the simulation as the systematic uncertainty in the background prediction as a function of p miss T and N b . As seen in Fig. 4, this uncertainty ranges from 8-15% in the lowest p miss T bin to 59-75% in the highest p miss T bin, and is assumed to be uncorrelated among bins.
Each of the listed uncertainties is incorporated in the background fit as a log-normal constraint in the likelihood function as described in Section 6.2, taking into account the stated correlations. Because of the robustness of the background prediction method, evidenced by the highstatistics data control region studies integrated in p miss T , the final uncertainty is dominated by the statistical precision of the simulation in evaluating the closure as a function of p miss T , described in the third item.

Results and interpretation
The observed event yields in data and the total predicted SM background are listed in Table 2, along with the expected yields for three higgsino mass scenarios. Two background estimates are given: the predictive fit, which does not use the data in the signal regions and ignores signal contamination in the other regions, and the global fit with r = 0, which incorporates the observations in the (HIG, 3b) and (HIG, 4b) regions, as described in Section 6.2. Since for p miss T > 450 GeV we observe no events in the (SBD, 4b) region, the parameter µ bkg 4b,SBD is fitted to be zero, pushing against its physical boundary and leading to the underestimation of the associated uncertainty. We account for this by including an additional contribution that makes the uncertainty on µ bkg 4b,SBD for p miss T > 450 GeV consistent with having observed one event. The event yields observed in data are consistent with the background predictions for all the analysis bins and no pattern of deviations is evident. The absence of excess event yields in data is interpreted in the context of the higgsino simplified model discussed in Section 1. Table 3 shows typical values for the systematic uncertainties associated with the expected signal yields for three models with different higgsino masses. The ranges correspond to the full variation of the uncertainties across all search bins. The uncertainty due to the pileup modeling is given by the difference between the signal efficiencies evaluated in samples with the mean number of reconstructed vertices found in the simulation and in the data, with the latter efficiencies obtained by extrapolation. The evaluation of the pileup uncertainty for very low higgsino masses is limited by the statistical power of the simulated samples. The remaining uncertainties are determined by comparing the nominal signal       yield for each search region to the corresponding yield obtained after varying the scale factor or correction under study within its uncertainty. In the case of the ISR uncertainty, the variation is based on the full size of the ISR correction derived by comparing the transverse momentum of the jet system balancing the Z boson in Z → + − events in data and in simulation. The largest uncertainties arise from the jet energy corrections, jet energy resolution, pileup modeling, and the p miss T resolution in the fast simulation. These uncertainties can be as large as 30% for low higgsino masses, but their impact is smaller for larger values of the higgsino mass. Uncertainties associated with the modeling of the b tagging range from 1% to 13%. The uncertainties in the trigger efficiency range from 6% in the lowest p miss T bin to <1% for p miss T > 300 GeV. Uncertainties due to the modeling of ISR and the efficiency of the jet identification filter are 1-2%. Finally, the systematic uncertainty in the total integrated luminosity is 2.5% [88].
The 95% confidence level (CL) upper limit on the production cross section for a pair of higgsinos in the context of the TChiHH simplified model is estimated using the modified frequentist CL S method [89][90][91], with a one-sided profile likelihood ratio test statistic in its asymptotic approximation [92]. Figure 9 shows the expected and observed exclusion limits. The theoretical cross section at NLO+NLL [45,46] as a function of higgsino mass is shown as a solid red line and the corresponding uncertainty as a dotted red line. The upper limits on the cross section at 95% CL for each mass point are obtained from the global fit method, which takes into account the expected signal contribution in all of the bins. Higgsinos with masses between 230 and 770 GeV are excluded.  Figure 9: Expected (dashed black line) and observed (solid black line) excluded cross sections at 95% CL as a function of the higgsino mass. The theoretical cross section for the TChiHH simplified model is shown as the red solid line.

Summary
A search for an excess of events is performed in proton-proton collisions in the channel with two Higgs bosons and large missing transverse momentum (p miss T ), with each of the Higgs bosons reconstructed in its H → bb decay. The data sample corresponds to an integrated luminosity of 35.9 fb −1 at √ s = 13 TeV. Because the signal has four b quarks, while the background is dominated by tt events containing only two b quarks from the t quark decays, the analysis is binned in the number of b-tagged jets. In each event, the mass difference between the two Higgs boson candidates is required to be small, and the average mass of the two candidates is used in conjunction with the number of observed b tags to define signal and sideband regions. The observed event yields in these regions are used to obtain estimates for the standard model background in the signal regions without input from simulated event samples. The data are also binned in regions of p miss T to enhance the sensitivity to the signal.
The observed event yields in the signal regions are consistent with the background predictions. These results are interpreted in the context of a model in which each higgsino decays into a Higgs boson and a nearly massless lightest supersymmetric particle (LSP), which is weakly interacting. Such a scenario occurs in gauge-mediated supersymmetry breaking models, in which the LSP is a goldstino. The cross section calculation assumes that the higgsino sector is mass degenerate and sums over the cross sections for the pair production of all relevant combinations of higgsinos, but all decays are assumed to be prompt. Higgsinos with masses in the range 230 to 770 GeV are excluded at 95% confidence level. These results constitute the most stringent exclusion limits on this model to date.
CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Aus   [93] CMS Collaboration, "Search for supersymmetry with Higgs boson to diphoton decays using the razor variables at √ s = 13 TeV", (2017). arXiv:1709.00384. Submitted to Phys. Lett. B.