Search for higgsinos decaying to two Higgs bosons and missing transverse momentum in proton-proton collisions at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 13 TeV

Results are presented from a search for physics beyond the standard model in proton-proton collisions at s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s} $$\end{document} = 13 TeV in channels with two Higgs bosons, each decaying via the process H → bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{b}} $$\end{document}, and large missing transverse momentum. The search uses a data sample corresponding to an integrated luminosity of 137 fb−1 collected by the CMS experiment at the CERN LHC. The search is motivated by models of supersymmetry that predict the production of neutralinos, the neutral partners of the electroweak gauge and Higgs bosons. The observed event yields in the signal regions are found to be consistent with the standard model background expectations. The results are interpreted using simplified models of supersymmetry. For the electroweak production of nearly mass-degenerate higgsinos, each of whose decay chains yields a neutralino χ~10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left({\overset{\sim }{\upchi}}_1^0\right) $$\end{document} that in turn decays to a massless goldstino and a Higgs boson, χ~10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left({\overset{\sim }{\upchi}}_1^0\right) $$\end{document} masses in the range 175 to 1025 GeV are excluded at 95% confidence level. For the strong production of gluino pairs decaying via a slightly lighter χ~20\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left({\overset{\sim }{\upchi}}_2^0\right) $$\end{document} to H and a light χ~10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left({\overset{\sim }{\upchi}}_1^0\right) $$\end{document}, gluino masses below 2330 GeV are excluded.


Introduction
The discovery of the Higgs boson (H) at the CERN Large Hadron Collider (LHC) [1][2][3][4] provides a new tool for probing physics beyond the standard model (SM). Higgs boson production is expected to be prominent in the decays of a variety of new, heavy particles, such as those arising in theories based on supersymmetry (SUSY) [5][6][7][8][9][10][11][12][13][14]. For example, Higgs bosons would be produced in the decays of higgsinos, their supersymmetric partners. A key motivation to search for higgsinos is that in certain models they are generically expected to be among the lightest supersymmetric particles, with masses kinematically accessible at the LHC. This is particularly the case in so-called natural SUSY models [15][16][17], which are designed to stabilize the electroweak mass scale. Searches for higgsinos are nevertheless challenging, not only because of the wide range of possible physics scenarios, but also because their electroweak coupling to other particles leads to small production cross sections. Higgsino searches are therefore expected to be a long-term part of the LHC physics program.

JHEP05(2022)014
We search for processes in proton-proton (pp) collisions leading to the production of two Higgs bosons, together with missing momentum transverse to the beam direction ( p miss T ) and possible additional jets. The Higgs bosons are reconstructed in the decay H → bb, which has a branching fraction of 58% [18]. To provide sensitivity to the wide range of kinematic configurations that can arise, the analysis searches both for events containing pairs of H → bb decays in which the b-tagged jets are separately resolved, and for events in which the b jets from a given H boson decay are merged together into a single, wider jet. We refer to the unmerged and merged kinematic configurations as the resolved signature and the boosted signature, respectively, because the merging of b jets typically occurs when the parent Higgs boson has a high momentum.
[51] for the boosted signature. Tabulated results are provided in the HEPData record for this analysis [65].
The simplest SUSY framework is the minimal supersymmetric standard model (MSSM), which includes an additional doublet of complex scalar Higgs fields [12][13][14] beyond the doublet in the SM. The SUSY partners of the gauge and Higgs bosons, referred to as gauginos and higgsinos, respectively, are all spin J = 1/2 fermions. In the MSSM, there are four higgsinos, two of which are charged ( H ± ) and two of which are neutral ( h 0 , H 0 ). However, the phenomenology of the electroweak sector of the MSSM is complicated by the fact that the partners of the electroweak gauge and Higgs bosons can mix. The resulting physical states are then the electrically neutral mass eigenstates, χ 0 1-4 (where the particles are ordered in mass, with χ 0 1 corresponding to the lightest), referred to as neutralinos, and the electrically charged mass eigenstates, χ ± 1,2 [66], referred to as charginos. In models that conserve R-parity [11,67], the lightest supersymmetric particle (LSP) is stable, and because of its weak interactions would escape experimental detection. The χ 0 1 would be a candidate for the LSP, as would the goldstino ( G), a Nambu-Goldstone particle associated with the spontaneous breaking of global supersymmetry. Figure 1 shows three models that can lead to the experimental signature considered in this analysis. These scenarios are described using the framework of simplified SUSY models [68][69][70][71], in which many of the SUSY particles are assumed to be decoupled or otherwise irrelevant to the process under consideration. In the TChiHH-G simplified model ( figure 1, left), the LSP is the goldstino. In a broad range of scenarios in which SUSY breaking is mediated at a low scale, such as gauge-mediated supersymmetry breaking (GMSB) models [72,73], the goldstino is nearly massless on the scale of the other particles and is the LSP. The χ 0 1 is then the next-to-lightest supersymmetric particle (NLSP) [74]. The NLSPs are produced in the cascade decays of several different combinations of neutralinos and charginos, and the goldstino is taken to be approximately massless. An important case arises when the lighter neutralinos χ 0 1,2 and charginos χ ± 1 are dominated by their higgsino content and, as a consequence, are nearly mass degenerate. In this case, all of their cascade decays can lead to the production of the NLSP and soft particles. The NLSP subsequently NLSPs are produced indirectly through the cascade decays of several combinations of neutralinos and charginos, as described in the text; (center) TChiHH, in which the electroweak production of two neutralinos leads to two Higgs bosons and two neutralinos ( χ 0 1 ); (right) T5HH, the strong production of a pair of gluinos, each of which decays via a three-body process to quarks and a neutralino, with the neutralino subsequently decaying to a Higgs boson and a χ 0 1 LSP. In each diagram, the hatched circle represents the sum of processes that can lead to the SUSY particles shown. decays to the goldstino and a Higgs boson, with a coupling that is suppressed by the SUSY breaking scale [72]. Integrating over the contributions from the allowed combinations of produced charginos and neutralinos ( χ 0 ) leads to an effective rate for χ 0 1 χ 0 1 production [75,76] that is significantly larger than that for any of the individual primary pairs. We assume a branching fraction of 100% for χ 0 1 → H G. In the TChiHH simplified model (figure 1, center), two higgsinos χ 0 2 and χ 0 3 are produced. The χ 0 1 is the LSP, assumed to be a bino (the superpartner of the SM boson corresponding to the U(1) weak hypercharge gauge field B), while χ 0 2 and χ 0 3 , nearly degenerate in mass, are the NLSPs. The other higgsinos have allowed decay channels, such as χ ± 1 → W ± χ 0 1 , that do not lead to Higgs bosons. This type of mass hierarchy can arise in various scenarios, as discussed in refs. [77,78]. Unlike in the TChiHH-G simplified model described above, where the heavier higgsino states were assumed to decay directly to the lightest one, here only the exclusive χ 0 2 χ 0 3 higgsino production cross section contributes; we are not sensitive to the other higgsino production channels χ ± 1 χ are not produced because their couplings to the Z boson are suppressed [78].) The cross section for χ 0 2 χ 0 3 alone is about 17% of that for the sum of all six higgsino cross sections [75,76].
Some SUSY models [79,80] predict production rates for energetic Higgs bosons that are greater in gluino cascade decays than in direct higgsino production, because of the much larger strong production cross section for gluinos. Figure 1 (right) corresponds to a model (T5HH) in which two gluinos are produced, each of which decays via a three-body process to a quark, an antiquark, and a χ 0 2 . The χ 0 2 decays to a Higgs boson and a χ 0 1 , which is taken to be the LSP. The same signature also arises in a recently proposed model [81,82], motivated by dark-matter considerations, where a bino-like χ 0 3 becomes the NLSP, decaying via a Higgs boson to the χ 0 1 LSP. This paper is organized as follows. Section 2 provides an overview of the analysis strategy, while sections 3 and 4 describe the CMS detector and the simulated event samples, respectively. The event triggers and reconstruction of the data are presented in section 5, while event selection and the reconstruction of Higgs boson candidates are discussed -3 -

JHEP05(2022)014
in section 6. The background estimation is described in section 7, and the systematic uncertainties are discussed in section 8. Section 9 presents the results and interpretation, and section 10 gives a summary.

JHEP05(2022)014
(pileup) occurring in the same bunch crossing, or in adjacent crossings within the time resolution of the data acquisition. The number of collisions varies with the instantaneous luminosity, averaging about 29 over the data set [86].

Simulated event samples
While the evaluation of the SM background, described in section 7, is based primarily on observed event yields in multiple control regions in data, Monte Carlo (MC) simulated event samples nevertheless play an important role in the analysis. Such samples are used to optimize the event selection requirements and background prediction methods without biasing the result; to evaluate correction factors, typically near unity, to the background predictions; and to evaluate the acceptance and efficiency for signal processes.
We use the MadGraph5_amc@nlo 2.2.2 (2.4.2) [87,88] event generator with leading order (LO) precision for simulation of the SM production of tt, W+jets, Z+jets, and QCD processes, for the 2016 (2017-2018) data taking periods. The tt events are generated with up to three additional partons in the matrix element calculations. The W+jets and Z+jets events are generated with up to four additional partons. For the MadGraph LO samples, MLM parton matching is used [88]. Single top quark events produced through the s channel, diboson events such as those originating from WW, ZZ, or ZH production, and rare events such as those from ttW, ttZ, and WWZ production, are generated with MadGraph5_amc@nlo at next-to-leading order (NLO), with FXFX matching [89], except that WW events in which both W bosons decay leptonically are generated using the powheg v2.0 [90][91][92][93][94] program at NLO. This same powheg generator is used to describe both tW production and t-channel production of single top quark events at NLO precision.
Samples of simulated signal events are generated at LO using MadGraph5_amc@nlo 2.4.2, with up to two additional partons included in the matrix element calculations. The production cross sections for the electroweak models are computed at NLO plus nextto-leading-log (NLL) accuracy [75,76], while those for strong production are determined with approximate next-to-NLO (NNLO) plus next-to-next-to-leading logarithmic (NNLL) accuracy [95][96][97][98][99][100][101][102][103][104][105][106]. For the TChiHH-G model, these samples are generated by performing a scan of m( χ 0 1 ), with m( G) = 1 GeV. The TChiHH model is represented by a two-dimensional scan of m( χ 0 1 ) and m( χ 0 2,3 ). The TChiHH-G and TChiHH simulated samples are generated with both Higgs bosons constrained to decay via H → bb, and each event is weighted by the square of the standard model branching fraction. Auxiliary simulations show that inclusion of the non-bb decay modes of the Higgs boson would slightly improve the sensitivity, mainly for the boosted signature, but these additional modes are not included here. Additional details on the simulation of electroweak production models are given in ref. [45]. Events with gluino pair production (T5HH) are generated for a range of m( g ) values. The gluino decay is simulated with a three-body phase space model [107]. The mass of the intermediate χ 0 2 is set 50 GeV below m( g ), ensuring that the daughter Higgs bosons have a large Lorentz boost, m( χ 0 1 ) is set to 1 GeV, and all decay modes of the Higgs bosons are included. All simulated samples make use of the pythia 8.205 [108] program to describe parton showering and hadronization. The CUETP8M1 [109] (CP5 [110]) pythia tune was used -6 -
The detector response is modeled with Geant4 [113]. Normalization of the simulated background samples is performed using the most accurate cross section calculations available [87,93,94,[114][115][116][117][118][119][120][121][122], which generally correspond to NLO or NNLO precision. Because scans over numerous mass points are required for the signal models, the detector response for such events (except for those of the T5HH model) is described using the CMS fast simulation program [107,123], which yields results that are generally consistent with those from the simulation based on Geant4. Pileup pp interactions are superimposed on the generated events, with a number distribution that is adjusted to match the pileup distribution measured in data.

Triggers and event reconstruction
The data sample for the analysis region was obtained using triggers [84,85] that require p miss T, trig to exceed a threshold that varied between 90 and 140 GeV, depending on the LHC instantaneous luminosity. Single-lepton control samples, defined below, are collected with single-electron, single-muon, and p miss T, trig triggers, while dilepton control samples are collected with single-electron and single-muon triggers. The value of p miss T, trig is computed with triggerlevel quantities and therefore has somewhat poorer resolution than p miss T . The trigger efficiency for the analysis sample is measured as a function of p miss T , and of H T , the scalar sum of the transverse momentum (p T ) of jets, separately for each year of data taking. For the single-lepton control sample, it is measured as a function of p miss T , H T , and lepton p T , and for the dilepton control sample, as a function of the leading-lepton p T . For the analysis sample, the efficiency at p miss T = 150 GeV is 30-70%, depending on the trigger thresholds (which varied with the instantaneous luminosity); the efficiency rises with p miss T and is over 99% for p miss T > 260 GeV. The analysis of the boosted signature requires p miss T > 300 GeV and thus operates on the efficiency plateau; the resolved-signature analysis operates partly on the turn-on section of the efficiency curve, with a lower p miss T threshold of 150 GeV. These measured efficiencies are applied as weights to the simulated SM and signal events, for the analysis region and the control samples.
The reconstruction of physics objects in an event proceeds from the candidate particles identified by the particle-flow (PF) algorithm [124], which uses information from the tracker, calorimeters, and muon systems to identify the candidates as charged or neutral hadrons, photons, electrons, or muons.
The primary pp interaction vertex is taken to be the reconstructed vertex with the largest value of summed physics-object p 2 T , as described in section 9.4.1 of ref. [125]. The physics objects for this purpose are the jets (reconstructed using the anti-k T jet finding algorithm [126,127] with the charged particle tracks assigned to the vertex), isolated tracks (including those identified as leptons), and the associated missing transverse momentum -7 -

JHEP05(2022)014
(computed as the negative vector sum of the p T values of those objects). Charged particle tracks associated with vertices other than the primary vertex are removed from further consideration. The primary vertex is required to lie within 24 cm of the center of the detector in the direction along the beam axis and within 2 cm in the plane transverse to that axis.
Jets from the primary pp interaction are formed from the charged PF candidates associated with the primary vertex, together with the neutral PF candidates, again using the anti-k T algorithm [126], as implemented in the FastJet package [127]. The clustering is performed with distance parameter R = 0.4 (AK4 jets) when optimized for jets containing the fragmentation products of a single parton, and with distance parameter R = 0.8 (AK8 jets) for jets containing multiple partons. No removal of overlap jets between the AK4 and AK8 collections is applied. Jet quality criteria [128,129] are imposed to select quark and gluon jets while rejecting those from spurious sources, such as electronics noise and detector malfunctions. The jet energies are corrected for the nonlinear response of the detector [130]. The estimated contribution to the jet p T of neutral PF candidates from pileup is removed with a correction based on the area of the jet and the average energy density of the event [131], for AK4 jets, and with the PUPPI technique [132], for AK8 jets. Jets are required to have |η| < 2.4; AK4 and AK8 jets are required to have p T > 30 GeV and p T > 300 GeV, respectively. Finally, AK4 jets that have PF constituents matched to an isolated electron or muon are removed. To improve the consistency of the fast simulation description with respect to that based on Geant4, we apply a small correction, about 1%, to account for differences in the efficiency of the jet quality requirements [128,129].
To improve the modeling of jets associated with initial-state radiation (ISR), the prediction of MadGraph5_amc@nlo is compared with data in a control sample enriched in tt events by the requirement of two light charged leptons (ee, µµ, or eµ) and two b-tagged jets. The number of all remaining jets in the event is denoted N ISR jet . A correction factor is applied to simulated tt and strong-production signal events so that the N ISR jet distribution agrees with that in data. The central value of the correction ranges from 0.92 for N ISR jet = 1 to 0.51 for N ISR jet ≥ 6. The correction is found not to be necessary for tt samples that are generated with the CP5 tune. For the electroweak signal samples, we account for the effect of ISR on the p T distribution of the system of SUSY particles, p ISR T , by reweighting the distribution based on studies of the transverse momentum of Z+jets events [133]. The reweighting factors range between 1.18 at p ISR T ∼ 125 GeV and 0.78 for p ISR T > 600 GeV. The identification of jets originating from b quarks (b tagging) is integral to the reconstruction of Higgs bosons in the analysis and is performed separately for AK4 and AK8 jets. For AK4 jets we use a version of the combined secondary vertex algorithm based on deep neural networks (DeepCSV [134]), making use of Loose (L), Medium (M), and Tight (T) working points, as discussed in section 6. The detection efficiencies for true b quark jets, as well as the mistag rates for jets originating from up, down, or strange quarks or from gluons (light-parton jets), and for charm-quark jets, are measured with data samples enriched in, respectively, tt, QCD, and a combination of W+c and semileptonic tt events [134]. Typical b jet efficiencies for the three DeepCSV working points are 85, 70, and 50%, with corresponding light-parton jet mistag rates of 10, 1, and 0.1%, and charm jet mistag rates of 40, 12, and 2.5%, for the L, M, and T working points, respectively. The The identification of AK8 jets containing two b quarks, used in the analysis of the boosted signature, is performed with a deep-learning-based double-b tagging algorithm (DeepDoubleBvL, mass-decorrelated [135]) at its loose working point. The efficiencies for tagging AK8 jets from H → bb and for those that contain only light partons are 90 and 5%, respectively [135]. The efficiency of the double-b tagging is checked with a data sample enriched in g → bb; the mistag rate is checked with data in a lower-p miss T sideband (200 < p miss T < 300 GeV) of the baseline selection for the boosted signature. For both single and double b tagging, the simulation is reweighted to compensate for differences with respect to data based on measurements of the b tagging efficiency and mistag rate for each working point. Additional corrections are applied for differences in the flavor tagging efficiencies between the fast simulation and Geant4, where applicable.
The reconstructed p miss T is computed as the negative vector sum of the p T of all PF candidates, adjusted for known detector effects by taking into account the jet energy corrections described in ref. [136]. Filters are applied to reject events with well-defined anomalous sources of p miss T , such as calorimeter noise, nonfunctioning calorimeter cells, beam particles externally scattered into the detector, and other effects [136]. For data generated with the fast simulation, corrections of 0-12% are applied to account for differences in the modeling of p miss T . Because the targeted signature is fully hadronic, we veto events with isolated charged lepton candidates as part of the definition of the analysis region. Conversely, we require isolated electrons and muons in the selection of certain control samples. Electrons are reconstructed from showers in the ECAL matched to tracks in the silicon detectors, taking into account the effects of bremsstrahlung in the material of the tracker. Additional criteria are imposed on the shower shape in the ECAL and on the ratio of energies associated with the electron candidate in the HCAL and ECAL [137]. We use the "cut-based veto" working point of the electron identification algorithm [137]. Muons are constructed from muon detector track segments matched to silicon detector tracks; we require candidates to satisfy the "medium" working point of the muon identification [138]. For both electrons and muons, the track is required to be consistent with the particle's origination at the primary vertex. Electron and muon candidates must also satisfy |η| < 2.5 and 2.4, respectively. Requirements on the lepton p T are imposed that differ between veto leptons, where efficiency is prioritized, and control-sample leptons, where higher purity is desired; we require p T > 10 GeV for veto leptons, and p T > 30 GeV for those used to select events in the single-lepton control samples. For the dilepton control samples, the threshold p T is 30 (20) GeV for the leading (subleading) lepton.
To select electrons and muons that originate in the decay of W and Z bosons, while suppressing those from the decays of hadrons and from particles misidentified as e or µ, both are required to be isolated from other PF candidates. The isolation criterion is based on the variable I, which is the scalar p T sum of charged-hadron, neutral-hadron, and photon PF candidates within a cone of radius ∆R ≡ (∆φ) 2 + (∆η) 2 around the lepton direction, after pileup subtraction, divided by the lepton p T , where φ is the azimuthal angle. The radius of the cone depends on the lepton p T : ∆R = 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 [139]. The isolation requirement is I < 0.1 (0.2) for electrons (muons).
The efficiency of electron and muon identification is measured with a tag-and-probe method applied to Z → + − samples in the data [137,138]. Yields of these leptons in simulation are corrected to match the observations in data.
The efficiency of the charged lepton veto can be compromised if the charged lepton is lost, either because it is not reconstructed or fails the lepton selection requirements (including the isolation and p T threshold requirements), or because it is a tau lepton that decays hadronically. To recover some of this veto efficiency, we also reject events with any additional isolated tracks corresponding to leptonic or hadronic PF candidates. To reduce the influence of tracks resulting from pileup, such isolated tracks are considered only if their distance of closest approach along the beam axis to a reconstructed vertex is smaller for the primary vertex than for any other vertex. The tracks are required to satisfy |η| < 2.5, as well as p T > 10 GeV, or >5 GeV for a PF electron or muon. The isolation requirement is that the scalar p T sum of all other charged particle tracks within a cone of radius 0.3 around the track direction, divided by the track p T , must be less than 0.2 if the track is identified as a PF electron or muon, less than 0.1 otherwise. To limit the veto tracks to those likely to originate in W boson leptonic decays, we also restrict the transverse mass, Here p T is the transverse momentum of the particle and ∆φ is the azimuthal angle between the particle momentum and p miss T .

Event selection and reconstruction of Higgs boson candidates
In the first step of the analysis procedure, as outlined in section 2, baseline selection requirements are applied to define the analysis regions for both the resolved and boosted signatures. These baseline selection requirements are described below. Three of the baseline selection requirements are based on quantities that are in common to both signatures and are defined as follows:

A minimum requirement on p miss
T is applied: p miss T > 150 GeV for the resolved signature and p miss T > 300 GeV for the boosted signature.
2. For both signatures, events are excluded if any veto leptons or isolated charged particle tracks are present. These veto objects are defined in section 5.
3. To suppress QCD background, events in either signature are rejected if any of the four highest p T AK4 jets is approximately aligned with the missing momentum vector p miss T . Such alignment is an indication that the observed p miss T in the event is a consequence of jet mismeasurement and is not genuine. This set of requirements, which we refer to collectively as the ∆φ requirement, can only be imposed in the plane transverse to the beam axis. An event is vetoed if any of the azimuthal separations is small, specifically, if ∆φ p miss The remaining baseline selection requirements depend on the individual details of the Higgs boson mass reconstruction and the counting of b jets, and we discuss these for the two signatures in turn.

JHEP05(2022)014
For the resolved signature, the following additional criteria are used to define the baseline requirements: 1. To control combinatorics in the Higgs boson reconstruction, events must contain either 4 or 5 AK4 jets. Besides the jets associated with the Higgs boson candidates, this requirement allows for one additional jet in the event, for example, from initial-state radiation.
2. Higgs boson reconstruction is performed as follows. The four jets with the highest b-tag discriminator values are selected and used to form two Higgs boson candidates, each one decaying into two b jets. There are three pairings among these four jets that could form the two Higgs boson candidates. To select one of them without biasing the mass distributions towards the true Higgs boson mass, we use the pair with the smallest absolute value of the mass difference, ∆m bb , between the two candidate masses The average mass m bb of these two Higgs boson candidates is required to satisfy the loose baseline requirement m bb < 200 GeV. This approach exploits the fact that signal events contain two particles of identical mass, without using the value of the Higgs boson mass itself. As such, it prevents the sculpting of an artificial peak in the background at m(H), so that a peak observed at that mass would be meaningful, and the combinatorial background in the m bb spectrum can be investigated in the subsequent analysis.
3. To further suppress background from tt lost-lepton events, we compute the angular separation, ∆R, between the b-tagged jets within each Higgs boson candidate. In tt lost-lepton events, one of the t quark decays produces a W boson that decays leptonically. The b jet from that t quark decay is then usually paired with a mistagged jet from the other t quark decay in the event to form a Higgs boson candidate. The value of ∆R for that candidate will then typically be large. This background is therefore reduced by requiring ∆R max < 2.2, where ∆R max is the larger of the ∆R values associated with the two H boson candidates.
4. We define a b-counting variable, N b , which is the number of b-tagged jets but using optimized selections based on tight, medium, and loose categories (see section 5) as follows: We refer to these categories as 2b, 3b, and 4b, but it should be remembered that the b tagging thresholds on individual jets become looser as N b increases. The baseline selection requires N b ≥ 2. the background estimation, while additional binning is performed in the quantities p miss T and ∆R max to improve the sensitivity of the analysis. In each plot, the baseline selection is applied to all variables other than the one shown, and, for all distributions except N b , an additional requirement N b ≥ 3 is used to select the signal-like parts of the overall analysis region. To facilitate comparison of the distributions between data and simulation, the distributions from simulation are scaled to match the area of the data, using factors listed in the caption. The representative signals show the expected peaking in m bb near the H boson mass, as well as toward lower values of ∆R max . The binning indicated by the vertical dotted lines is discussed in section 7. While these plots provide an overall picture of the background shapes and background composition in the analysis region, they do not give an accurate picture of the sensitivity of the analysis to the signal, because they do not yet reflect the final binning used to extract the signal.
The additional baseline criteria for the analysis of the boosted signature are as follows: 1. At least two AK8 jets, each with p T > 300 GeV. The p T threshold results in a high probability that all daughter particles from the Higgs boson decay are contained within the jet. The kinematic merging of the Higgs boson decay products reduces the combinatorial challenge encountered in the resolved signature, and so the number of additional AK8 jets is unrestricted, as is the number of AK4 jets. Such additional jets can be produced in the signal processes involving strongly interacting particles, such as the T5HH model described in section 1. They can also arise from initial-or final-state radiation.
2. The mass m J attributed to an AK8 jet is computed by the "soft drop" algorithm [140,141], in which soft wide-angle radiation is recursively removed from a jet. For signal events, m J is expected to show a peak near m(H). For the baseline selection, the two highest p T AK8 jets are required to satisfy the loose requirement 60 < m J < 260 GeV.
3. We define the variable N H as the number of AK8 jets passing the above criteria that are double-b tagged, that is, for which the value of the double-b tagging discriminator D bb exceeds the loose working point threshold of 0.7. The value of N H , 0, 1, or 2, is unrestricted in the baseline selection; it is used for the event classification discussed in section 7. Figure 3 shows the distributions for data and simulated event samples with the boosted signature in p miss T , and in p T , D bb , and m J for the two leading AK8 jets. The baseline requirements for the boosted signature are applied, except for those on the variable shown. The simulated SM background yields, scaled by the factor noted in the caption, are seen to describe the shape of the data well. Also shown are examples of SUSY models, which peak around the Higgs boson mass in m J , and at high values of D bb , for both jets.
In a small region of parameter space a candidate may be selected by both signatures. In this case it is assigned to the resolved selection. Slightly higher expected sensitivity is achieved with this choice than with the alternative. The efficiency of the combined resolved and boosted signatures for the targeted signal models ranges from 1 to 15%, depending on the values of the mass parameters in the models.

Background estimation
After applying the baseline event selections for the resolved and boosted signatures described in section 6, we define the analysis regions shown in figure 4. The background in a given signal region (SR) is estimated using the event yields observed in three types of neighboring regions: the sideband regions (SB), the control sideband regions (CSB), and the control signal region (CSR). For both signatures, these regions are constructed using two types of discriminating variables: 1. Variable(s) characterizing the masses of the H boson candidates. The resolved signature uses a single variable, m bb , while the boosted signature uses both m J1 and m J2 , where J 1 and J 2 are the p T -leading and subleading AK8 jets, respectively.

2.
A b-counting variable derived from the jet flavor tagging. For the resolved signature, the variable is N b , while for the boosted signature, it is N H .
For the resolved signature, signal events primarily have four b-tagged jets, but they can also populate the region with three b-tagged jets, so the signal regions are defined by requiring m bb to be in the Higgs boson mass window and either N b = 4 or N b = 3. For the boosted signature, signal events primarily have two identified double-b-tagged jets, but they can also populate the region with only one. The signal regions are therefore defined by requiring both m J1 and m J2 to be in the Higgs boson mass window and either N H = 1 or N H = 2.
The background event yield in a signal region is then estimated with the "ABCD method," which uses the measured event yields in three background-dominated regions to JHEP05(2022)014 predict the background in a signal region. With the regions defined as shown in figure 4, the estimated background yield is where N SB is the observed event yield in the H boson mass sideband for events satisfying the same b-counting criteria as those applied in the SR, while N CSR and N CSB are, respectively, the event yields in regions corresponding to signal and sideband regions with respect to the Higgs candidate mass variable, but in the background regions with respect to the b-counting criterion. The coefficient κ is a correction factor, typically near unity, that takes into account a potential correlation between the two types of discriminating variables. It is taken from simulation but is validated with separate control regions in data, as discussed in section 8.
For the resolved signature, figure 4 (left), the variable N b was defined in section 6, and the H boson mass window for the SR or CSR is given by 100 < m bb < 140 GeV, with the events both above and below this window constituting the SB or CSB. The SRs and SBs reside in b tagging regions N b = 3 and N b = 4, while the N b = 2 region contains the CSR and CSB. To increase the sensitivity to signal, these regions are further sorted into bins in two discriminating variables, ∆R max (0-1.1 and 1.1-2.2) and p miss T (150-200, 200-300, 300-400, and >400 GeV). The background determination by the ABCD method is performed separately for each of these eight bins, and for both the N b = 3 and N b = 4 SRs, for a total of 16 bins.
The corresponding ABCD regions for the boosted signature, figure 4 (right), are based on the number N H of double-b-tagged jets (defined in section 6), with events sorted into regions 0H, 1H, and 2H having N H = 0, 1, and 2, respectively. The SRs and SBs lie within the 1H and 2H regions, and CSR and CSB in the 0H region. The events in the SR and CSR have 95 < m J < 145 GeV for both AK8 jets, while those in the SB and CSB have at least one of the jets lying outside that mass window, as illustrated in the right-hand plot of figure 4. The background yield N pred SR, tot is determined for each of the 1H and 2H SRs by the ABCD approach, eq. (7.1). The events in these SRs are further sorted into three p miss T bins (300-500, 500-700, and >700 GeV), for a total of six bins. The distribution of the predicted background yield among these p miss T bins is discussed below. The factors κ in eq. (7.1) can differ from unity if the H candidate mass variable is correlated with the b-counting variable, such that the N SR /N SB and N CSR /N CSB yield ratios differ. We thus evaluate κ as a double ratio taking the yields in this case from simulation. A strength of this method is that many systematic uncertainties cancel in this double ratio. Figure 5 shows the values obtained for κ in each of the 16 analysis bins of the resolved signature. These corrections, with their uncertainties, are applied to the background yield predictions, as indicated in eq. (7.1). The determination of the values of κ for the boosted signature is illustrated in figure 6 (left), which shows the ratios from simulation N CSR /N CSB for the 0H, and N SR /N SB for the 1H  Figure 5. The double ratio κ = (N SR /N SB )/(N CSR /N CSB ) from the SM simulation for the N b = 3 and N b = 4 search samples for each (p miss T , ∆R max ) bin of the resolved signature. The value ∆ gives the deviation of κ from unity, and σ its relative statistical uncertainty. and 2H analysis regions. Dividing this ratio for the 1H (2H) region by the ratio for the 0H region we obtain the corresponding factor, with statistical uncertainties, κ = 1.02 ± 0.04 (0.94±0.17). As these are consistent with unity, we set κ = 1 in the calculation of the central values of the background yields, propagating the uncertainties in κ to these yields. For both signatures, systematic contributions to the κ uncertainties are evaluated as discussed in section 8.
For the boosted signature, the data yields are too small to permit a separate ABCD background calculation in each p miss T bin. Instead we combine the p miss T -integrated estimate N pred SR, tot for each of the 1H and 2H SRs with the binned p miss T distribution f pTmiss obtained from a separate data control region. This control region is a subset (0H+b) of the 0H region, comprising both the CSR and CSB regions but satisfying the further requirement of one tight b-tagged AK4 jet in addition to the two AK8 jets. The b tagging requirement serves to enrich this control region in tt events, so as to approximate the SM content of the 1H and 2H regions. We thus have for the predicted yield, N pred SR,i , in p miss where N pred SR, tot is the predicted background yield, as calculated from eq. (7.1) without any binning in p miss T , and f pTmiss,i is the fraction of events in the 0H+b region that fall in p miss T bin i. Specifically, and κ is a correction factor obtained from simulation, as described below. The yields N 0H+b,1,2,3 in the three p miss T bins of the 0H+b region are included in the likelihood function described in section 9 so as to propagate the statistical uncertainties of f pTmiss to the signal yield. We use the SM simulation to test the consistency of the p miss T distribution in the 0H+b region with those of the 1H and 2H regions. This test is shown in figure 6 (right). The ratio distributions in the lower panel are seen to be uniform within the uncertainties of the simulation. We thus set the central values of the κ factors to unity, accounting for their uncertainties as described in section 8. The figure also shows agreement between the simulation and the fractional distribution f pTmiss as measured in the data (open points).

Systematic uncertainties
Systematic uncertainties in the background estimates arise principally from the statistical uncertainties in control regions of the data, and from statistical and systematic uncertainties in the simulated samples used to measure κ and κ . These quantities can be measured in SM-dominated control samples, providing a check on the procedure and estimates of the associated systematic uncertainties. Uncertainties in the simulation also affect the signal predictions. Sections 8.1 and 8.2 describe the evaluation of the background systematic uncertainties for the resolved and boosted signatures, respectively. Section 8.3 describes the systematic uncertainties in the signal model predictions.

Background systematic uncertainties (resolved signature)
We validate the κ factors described in section 7 and displayed in figure 5, and assign systematic uncertainties, by comparing their values obtained in simulation with those from three control samples, each of which is enriched in a specific class of SM background process. We compare the value of κ computed from the data with the one computed from the corresponding simulation. Since for each control sample we find the comparison to be consistent across p miss T bins, we combine these bins, and take for the systematic uncertainty -18 -

JHEP05(2022)014
the larger of the data-simulation difference and the statistical uncertainty in the value obtained from simulation.
The tt systematic uncertainty is measured in a tt-dominated single-lepton control sample selected by replacing the lepton and charged particle vetoes of the search sample with the requirement of one electron or muon with p T > 30 GeV, and no additional muons or electrons. The ∆φ requirement is relaxed in the absence of QCD background, to improve the statistical precision, and a requirement m T < 100 GeV is imposed to suppress signal contamination. We measure the relative systematic uncertainties in κ for the tt background in the ∆R max < 1.1 bins to be 13% and 19% for N b = 3 and N b = 4, respectively; the corresponding numbers for the 1.1 < ∆R max < 2.2 bins are 2% and 8%.
The V+jets systematic uncertainty is measured in a Z+jets dominated dilepton control sample in which we select events with two electrons or muons of opposite sign, at least one of which has p T > 30 GeV, dilepton mass between 80 and 100 GeV, and p miss T < 50 GeV. We then compute the p T of the dilepton system, p T (Z), and treat this quantity as the true p miss T . That is, p T (Z) serves as a proxy for the p miss T that would be measured in an event with a Z → νν decay. To increase the sample size, we relax the ∆φ requirement and extend the range of p T (Z) down to zero. To reduce tt contamination, we take the CSR and CSB (SR and SB) regions to be defined by the b jet selection N b,M = 0 (1). We measure the relative systematic uncertainties in κ for the V+jets background to be 16% and 5% for the ∆R max < 1.1 and 1.1 < ∆R max < 2.2 bins, respectively.
The QCD multijet systematic uncertainty is measured in a QCD-enriched sample selected by inverting the ∆φ requirement of the search sample. As with the dilepton control sample, the CSR and CSB (SR and SB) are defined by the b jet selection N b,M = 0 (1). We measure the relative systematic uncertainties in κ for the QCD background to be 9% and 7% for the ∆R max < 1.1 and 1.1 < ∆R max < 2.2 bins, respectively.
Finally, these results for each SM process are weighted by the fraction of total background arising from that process in each of the 16 analysis bins to obtain the uncertainties given in table 1. This direct in situ comparison with data of the simulation affecting κ covers all of the background-related uncertainties associated with MC modeling.
The effect of trigger efficiency modeling on κ has been evaluated and found to be negligible (<0.1%). Table 1 gives a summary of the uncertainties in the κ factors for the resolved signature. These are propagated to the background prediction through eq. (7.1). The table shows that in most of the SRs the largest contribution to the systematic uncertainty in the background prediction is the statistical uncertainty in the simulation used to determine κ, which ranges from 4 to 56%.

Background systematic uncertainties (boosted signature)
Systematic uncertainties in the background yields are evaluated for each of the two steps (eq. (7.1) and eq. (7.3)) in the background determination for the boosted signature, and are summarized in  Table 1. Summary of the uncertainties in κ for the resolved signature. The sources are statistical uncertainties in the determination of κ and the systematic uncertainties ∆κ(data, MC)/κ derived from the comparison of simulation with data in the single-lepton, dilepton, and low-∆φ control samples, each weighted by the fraction of background arising from the associated process in each search bin.
by applying the baseline selection for the boosted-signature analysis region, but releasing the lepton and isolated track vetoes and the ∆φ requirement, and instead requiring exactly one electron or muon with p T > 30 GeV and m T < 100 GeV. Taking the larger of the datasimulation difference and the statistical uncertainty in that comparison, we find systematic uncertainties of 9 and 13% for the 1H and 2H SRs, respectively. The statistical uncertainty from the calculation of κ in simulation is assigned as a further systematic uncertainty. The second step in the background prediction for the boosted signature is the application of the p miss T distribution f pTmiss measured with the 0H+b control region. Systematic uncertainties arising from possible differences in the p miss T shape between the 0H+b control region and the SRs are contained in the correction factors κ in eq. (7.3). These uncertainties are measured as the deviations, from the p miss T -averaged value, of a linear fit to the data in the lower panel of figure 6 (right); the values are 2, 3, and 9% (13, 19, and 63%) for the three p miss T bins in the 1H (2H) SR. We obtain an additional uncertainty in the p miss T distribution, attributable to the background composition of the simulation, by varying the Z+jets component by factors of two and one half. Taking the larger effect of these on κ yields the uncertainties 3, 7, and 32%.
We assign a systematic uncertainty covering the effects of MC mismodeling in both N pred SR, tot and f pTmiss from the following sources: b and double-b tagging, jet energy scale, jet energy resolution, ISR corrections, and renormalization (µ R ) and factorization (µ F ) scales. The effect of each of these contributions ranges from 0 to 5%. The total systematic -20 -  Table 2. Summary of systematic uncertainties in the background prediction for the boosted signature. Values in cells spanning the 1H and 2H columns enter the yield calculations through the factor f pTmiss , which is common to the 1H and 2H SRs. The values from the ABCD measurement with the p T -integrated sample appear in the rows labeled p miss T > 300 GeV. The row labeled ∆(data, MC) gives the contribution derived from the comparison of simulation with data in the one-lepton control sample. uncertainty is calculated by taking the largest uncertainty for each source, and adding those in quadrature. The individual sources are described in section 5 and their uncertainties in section 8.3.

JHEP05(2022)014
As seen in table 2, the leading contributions are statistical uncertainties in the SB, CSR, CSB, and 0H+b yields in the data, and in those of the MC samples that are used to measure κ .

Signal systematic uncertainties
Systematic uncertainties affecting the signal yields from simulation are listed in table 3. The systematic uncertainty in the ISR correction of the simulation of electroweak production models is evaluated by varying the correction factors described in section 5 (18-22%) by 100% of their nominal value. For strong production models, the systematic uncertainty is evaluated by varying the correction factors (4-26%) by 50% of their nominal value. To evaluate the uncertainty associated with µ R and µ F , we vary one or both by factors of 2.0 and 0.5 and take the combination (excluding the case of opposite variations of the two) that  Table 3. Sources of systematic uncertainties and their typical impact on the signal yields obtained from simulation. The range is reported as the median 68% confidence interval among all signal regions for every signal mass point considered. Entries reported as 0 correspond to values smaller than 0.5%.
gives the largest effect on the yield [142][143][144]. The uncertainty associated with the pileup reweighting is evaluated by varying the value of the total inelastic cross section by 5% [145]. The systematic uncertainty in the determination of the integrated luminosity varies between 1.2 and 2.5% [146][147][148], depending on the year of data collection; the total integrated luminosity for the three years has an uncertainty of 1.6%. The uncertainties related to the jet energy scale and jet energy resolution are calculated by varying jet properties in bins of p T and η, according to the uncertainties in the jet energy corrections. The uncertainty in the m J resolution is derived by comparing the mean and width of the W boson peak in the m J spectrum between data and simulation in a tt-enriched sample. The efficiency of the isolated-track veto for signal events, as measured in simulation, is greater than 90% except at the lowest NLSP masses. Differences with respect to simulation measured in data control samples lie in the range 20-30% of the inefficiency. We assign a systematic uncertainty of 50% of the inefficiency to account for the effect of this veto on the signal yields.

JHEP05(2022)014
The uncertainty in the trigger efficiency modeling is derived from the statistical uncertainty, and variations in the kinematic selections, for the trigger measurement sample, and from variation in the independent trigger used to measure the efficiency. This uncertainty is applied to each SR bin for each signal model based on its distribution in p miss T and H T in that bin. For the resolved signature, typical values range from 1-13%, though some bins have larger values. For p miss T > 300 GeV the largest uncertainty is 4%. Systematic uncertainties in the b tagging are obtained from the measurements in data of the b and double-b tagging efficiencies and mistag rates, described in section 5. These uncertainties lie in the ranges 2-7% for jets tagged as b quarks and 6-15% for the double-b-tagged jets.
Systematic uncertainties associated with use of the fast simulation are taken to be 100% of the differences with respect to the Geant4 simulation. This amounts to 1% for the jet quality requirements, 2-4% for m J resolution, and 1% for b and double-b tagging.

Results and interpretation
The observed distributions in the appropriate Higgs candidate mass quantity ( m bb for the resolved signature and m J for the boosted signature) are shown in figures 7-9.  figure 9 shows the data distributions in m J , and in (m J1 , m J2 ), with simulation overlaid, for the 0H CSR+CSB and the 1H and 2H SR+SB regions, integrated in p miss T . Again, we do not observe a large excess in the data in the m J distributions at the Higgs-boson mass. Systematic uncertainties are not included in these figures, as those uncertainties are computed only for the total yields in the analysis regions. The plots also show representative distributions for potential SUSY signals, with their clear peaks in the H boson mass window.
To extract the signal strength, µ, we perform maximum likelihood fits to the data using a probability density function for each bin that represents the SM yield plus the yield predicted by a specific model for the signal. The parameter µ is the ratio of the fitted to theoretical cross section. The likelihood function is a product of Poisson distributions, one for each of the SRs of the resolved and boosted signatures and their corresponding SB, CSR, and CSB regions. Additional Poisson variables represent the event yields in the three p miss   Table 4. For each SR of the resolved signature, the MC correction factor κ, predicted background yield N pred SR , yield from the background-only (µ = 0) fit N fit SR , and observed yield N obs SR . The first and second uncertainties in the κ factors are statistical and systematic, respectively. The uncertainties in N pred SR and N fit SR , extracted from the maximum likelihood fit, include both statistical and systematic contributions. The interpretation of the results in bin 11 is discussed in the text. component of the Poisson means for each of these factors. The measured values of κ, with their statistical uncertainties, are introduced as Gaussian constraints. Other systematic uncertainties are implemented as additional free parameters with log-normal constraints. Correlations among signal bins are taken into account. Inclusion of the signal component in the functions for the Poisson means of the SB, CSR, CSB, and 0H+b regions accounts for potential signal contamination of these regions.
We determine confidence intervals for µ using the test statistic q µ = −2 ln(L µ /L max ), where L max is the maximum likelihood determined by allowing all parameters including µ to vary, and L µ is the maximum likelihood for a fixed signal strength. Confidence ranges are set under the asymptotic approximation [149], with q µ approximated with an Asimov data set [149] and used in conjunction with the CL s criterion described in refs. [150,151].
The observed yields, together with the predicted SM yields, are given in tables 4 and 5 for the resolved and boosted signatures, respectively, and are summarized in figure 10.  Table 5. For each N H SR of the boosted signature, the total predicted background yield N pred SR, tot , and for each p miss T bin the fraction f pTmiss , both with their statistical uncertainties, the predicted background yield N pred SR , the yield from the background-only (µ = 0) fit N fit SR , and the observed yield N obs SR . The values of N pred SR and N fit SR are extracted from the maximum likelihood fit, with uncertainties that include both statistical and systematic contributions.  the background-only fit (column N fit SR in table 4) the background estimate becomes 0.72 +0.53 −0.33 . We have examined the four observed events in detail and find no exceptional features. The models we consider for potential contributions beyond the SM predict that events would be distributed over multiple bins, contrary to the observations here. To account for the look-elsewhere effect [152], we perform pseudo-experiments with all 16 bins of the resolved signature. We find that the probability of observing a 3.3 s.d. or larger excess in at least one bin corresponds to a global significance of 2.1 s.d.
We interpret the data in terms of the three models discussed in section 1. In each case, the data from both resolved and boosted signatures are included in the fit. Figure 11 shows the upper limit on the cross section at 95% CL for the TChiHH-G model of figure 1 (left). For the case of the production of the sum of all higgsino pairs, the data exclude masses of the χ 0 1 in the range 175-1025 GeV. This represents the best limit to date on this model, extending previous ATLAS results [29], which exclude mass ranges 130-  For the TChiHH model (figure 1, center) the cross section limit is presented in figure 13 as a function of the independent masses m( χ 0 2 ) and m( χ 0 1 ). While the expected limit would exclude a substantial region of this plane, the observed exclusion is limited to a small region where m( χ 0 1 ) is less than 15 GeV, again because of the single-bin excess. For m( χ 0 1 ) = 1 GeV the excluded range of m( χ 0 2 ) is 265-305 GeV. Figure 14 shows the cross section upper limit as a function of m( g ) for the T5HH model ( figure 1, right) of gluino pair production with a χ 0 2 NLSP slightly less massive than the gluino and a light LSP. Masses m( g ) < 2330 GeV are excluded. This is the strongest mass limit for this model to date, extending a previous CMS result [51], which excluded m( g ) < 2010 GeV. Most of the sensitivity to this model is provided by the boosted signature. This reflects the choice of large NLSP mass, which leads to energetic H boson daughters.

Summary
A search has been presented for physics beyond the standard model in channels leading to pairs of Higgs bosons and an imbalance of transverse momentum, produced in proton-proton collisions at √ s = 13 TeV. The data sample, collected by the CMS experiment at the LHC, corresponds to an integrated luminosity of 137 fb −1 . The Higgs bosons are reconstructed via their decay to a pair of b quarks, observed either as distinct b quark jets (resolved signature), or as wide-angle jets that each contain the pair of b quarks (boosted signature). The observed event yields in 15 of the 16 analysis bins of the resolved signature, and all 6 bins of the boosted signature, are consistent with the background predictions based on SM processes. In one bin of the resolved signature, an excess is observed, for which the global significance is 2.1 standard deviations when all 16 bins are considered. A narrow, single-bin effect is not typical of the supersymmetry (SUSY) models under consideration.
These results are used to set limits on the cross sections for the production of particles predicted by SUSY, considering both the direct production of neutralinos and their production through intermediate states with gluinos. For the electroweak production of nearly-degenerate higgsinos, each of whose decay cascades yields a χ 0 1 that in turn decays Observed and expected upper limits at 95% CL on the cross section for the simplified model T5HH, the strong production of a pair of gluinos each of which decays via a three-body process to quarks and a χ 0 2 NLSP, which subsequently decays to a Higgs boson and a χ 0 1 LSP. The dashed black line with green and yellow bands shows the expected limit with its 1-and 2-s.d. uncertainties; the solid black line shows the observed limit, and the dashed red line the theoretical cross section. through a Higgs boson to the lightest SUSY particle (LSP), a massless goldstino, χ 0 1 masses in the range 175-1025 GeV are excluded at 95% confidence level. For a model with a mass splitting between the directly produced, degenerate higgsinos χ 0 2 χ 0 3 , and a bino LSP, a small region where m( χ 0 1 ) is less than 15 GeV is excluded; for m( χ 0 1 ) = 1 GeV the excluded range of m( χ 0 2 ) is 265-305 GeV. For the strong production of gluino pairs decaying via a slightly lighter χ 0 2 to a Higgs boson and a light χ 0 1 LSP, gluino masses below 2330 GeV are excluded. The bounds on masses found here extend previous limits for these models by about 150 and 320 GeV for the χ 0 1 and gluino, respectively.