Search for light Higgs bosons from supersymmetric cascade decays in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {pp}$$\end{document}pp collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13\,\textrm{TeV}$$\end{document}s=13TeV

A search is reported for pairs of light Higgs bosons (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textrm{H}} _1$$\end{document}H1) produced in supersymmetric cascade decays in final states with small missing transverse momentum. A data set of LHC \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {pp}$$\end{document}pp collisions collected with the CMS detector at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13\,\text {TeV} $$\end{document}s=13TeV and corresponding to an integrated luminosity of 138\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {fb}^{-1}$$\end{document}fb-1 is used. The search targets events where both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textrm{H}} _1$$\end{document}H1 bosons decay into pairs that are reconstructed as large-radius jets using substructure techniques. No evidence is found for an excess of events beyond the background expectations of the standard model (SM). Results from the search are interpreted in the next-to-minimal supersymmetric extension of the SM, where a “singlino” of small mass leads to squark and gluino cascade decays that can predominantly end in a highly Lorentz-boosted singlet-like \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textrm{H}} _1$$\end{document}H1 and a singlino-like neutralino of small transverse momentum. Upper limits are set on the product of the squark or gluino pair production cross section and the square of the branching fraction of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textrm{H}} _1$$\end{document}H1 in a benchmark model containing almost mass-degenerate gluinos and light-flavour squarks. Under the assumption of an SM-like branching fraction, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textrm{H}} _1$$\end{document}H1 bosons with masses in the range 40–120\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {GeV}$$\end{document}GeV arising from the decays of squarks or gluinos with a mass of 1200–2500\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\,\text {GeV}$$\end{document}GeV are excluded at 95% confidence level.


Introduction
This paper presents a search for pairs of light Higgs bosons (H 1 ) produced in supersymmetric (SUSY) [1][2][3][4][5][6][7][8] cascade decays in final states with small missing transverse momentum (p miss T ). Such events can arise from the pair production of squarks ( q) and gluinos ( g) in the next-tominimal supersymmetric extension of the standard model (SM) [9] when the lightest SUSY particle (LSP) is a singlino-like neutralino ( χ 0 S ) of small mass [10]. The χ 0 S mass eigenstate is dominated by the singlino component and has only small couplings to other SUSY particles, suppressing direct squark or gluino decays to the χ 0 S . Squarks and gluinos decay via the nextto-LSP χ 0 2 into a χ 0 S and a Higgs, Z, or W boson [10,11]. The case of a singlet-like CP-even H 1 , shown in Fig. 1, is the focus of this search. When the χ 0 S has a far smaller mass than the H 1 and the phase space for the decay χ 0 2 → H 1 + χ 0 S is small, the H 1 carries much larger momentum than the χ 0 S . In such p miss T -suppressed scenarios, the key signature for the pair production of squarks and gluinos is a pair of Lorentz-boosted H 1 bosons.  Figure 1: Diagram of squark pair production and subsequent cascade decay in the benchmark signal model. The particle χ 0 2 is the next-to-LSP, χ 0 S is the singlino-like LSP, and H 1 is the CPeven singlet-like Higgs boson.
This search targets events with two highly boosted H 1 bosons that decay into bb pairs that are reconstructed as large-radius jets using substructure techniques. This is the first search at the LHC to focus on this type of event, where particles invisible to the detector have only small transverse momentum (p T ) and therefore the events are not selected by searches requiring significant p miss T [10,12]

Benchmark signal model
A benchmark signal model is established following the work of Ellwanger and Teixeira [10,11]. The eight first-and second-generation squarks are assumed mass-degenerate at the mass m SUSY , and the gluino mass is set at 1% larger. The small gluino-squark mass gap means that the kinematics of the final-state particles are very similar in the q q, q g, and g g production modes, as little momentum is transferred to the quark in the g → q + q decay. All SUSY particles other than gluinos and those shown in Fig. 1 are assumed decoupled. This search targets squarks and gluinos with m SUSY > 1200 GeV. Less massive squarks and gluinos can be probed by p miss T -based searches, owing to their larger pair-production cross sections [12]. Smaller m SUSY values can also lead to smaller p T of the H 1 than is necessary for the bb pair to be merged in a single jet. The cross sections (σ) for the signal probed in this search, calculated at next-to-leading order (NLO) accuracy in the strong coupling constant (α S ) including approximate next-to-NLO (NNLO) corrections and next-to-next-to-leading logarithmic (NNLL) soft gluon corrections [21][22][23][24][25][26][27][28][29], are shown in Table 1. Table 1: Inclusive pair-production cross sections calculated at approximately NNLO and NNLL in α S [21][22][23][24][25][26][27][28][29] for squark mass m SUSY and gluino mass 1% larger. The quoted uncertainty is obtained from variations in the choice of scales, parton distribution functions, and α S . m SUSY [GeV] σ(pp → q q, q g, g g ) [ [10], as calculated using HDECAY 6.61 [30,31]. The B values decrease for larger H 1 masses as the virtual WW * and ZZ * decay channels, both of which have sizeable leptonic branching fractions, become more accessible. The region m H 1 < m Z is therefore where the p miss T -suppressed all-jet signature is of greatest experimental importance. Nevertheless, to preserve generality, this search attempts to probe as much of the region m H 1 < 125 GeV as possible.
The p miss T -suppressed signature arises for values of R m close to unity, provided ∆ m > 0 to permit the χ 0 2 → H 1 + χ 0 S decay. In this case, the phase space for the χ 0 2 decay is small and the χ 0 S has much smaller mass than the H 1 , so the χ 0 S always carries much less momentum than the H 1 . The p miss T -suppressed signature probed in this search is representative of a significant part of the model parameter space since the momenta of reconstructed objects do not exhibit a strong dependence on R m and ∆ m in the region R m > 0.9. Models with smaller R m can be probed by p miss T -based searches [10,12]. For the benchmark model, the values R m = 0.99 and ∆ m = 0.1 GeV are assumed.
Branching fractions of unity are assumed for the decays q → q + χ 0 2 and χ 0 2 → H 1 + χ 0 S . In the R m and ∆ m region of the benchmark model, this is true except where m In that case, the χ 0 2 → Z + χ 0 S decay is permitted if the χ 0 2 has a higgsino component [11]. However, the χ 0 2 is expected to be mainly bino-like for relevant values of its mass [10]. For configurations where the H 1 mass is close to that of the H SM , the decay χ 0 2 → H SM + χ 0 S is also possible. The signatures for such H 1 and H SM bosons are indistinguishable in this search. The assumption that the branching fraction for χ 0 2 → H 1 + χ 0 S decay is 100% can therefore be relaxed to the assumption that the branching fractions to H 1 and H SM sum to unity.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. A silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections, reside within the solenoid volume. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are measured in gas-ionisation detectors embedded in the steel flux-return yoke outside the solenoid. Events of interest are selected using a two-tiered trigger system. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed latency of about 4 µs [32]. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimised for fast processing, and reduces the event rate to around 1 kHz before data storage [33]. A more detailed description of the CMS detector, together with a definition of the coordinate system and the kinematic variables, can be found in Ref. [34].

Event simulation
The primary background in this search originates from multijet production. Simulated multijet events are used to validate the multijet background estimation based on data (described in Section 6), but are not used for any of the final predictions. The remaining significant background is from events with vector bosons that decay into quark-antiquark pairs. Simulated events are used to determine the contributions from tt, Z+jets, and W+jets production. The expected yields from all other SM sources of background are found to be negligible.
The multijet, Z+jets, and W+jets processes are simulated at leading order (LO) in perturbative quantum chromodynamics (QCD) using MADGRAPH5 aMC@NLO 2.4.2 [35] with up to four additional partons at the matrix element (ME) level. Simulated signal events for each pair of m H 1 and m SUSY values of the benchmark model are generated at LO at the ME level with up to one additional parton using MADGRAPH5 aMC@NLO 2.3.3. The MLM [36] prescription is used to match partons from the LO ME calculations to those from the parton showers. Simulated tt events are produced at NLO in QCD at the ME level with the POWHEG v2.0 [37][38][39][40] generator. The NNPDF2.3, NNPDF3.0, and NNPDF3.1 [41][42][43][44] parton distribution functions (PDFs) are used for the signal, 2016 background, and 2017-2018 background simulations, respectively. The parton shower and hadronisation are performed via PYTHIA 8.2 [45]. The CUETP8M1 [46,47] tune is used for the signal and 2016 background simulations, while the CP5 tune [48] is used for the 2017 and 2018 background simulations. The cross section used to normalise the tt simulation is calculated at NNLO+NNLL in QCD [49], and those for Z+jets and W+jets are calculated at NNLO in QCD [50][51][52]. Additional pp interactions within the same or nearby bunch crossings (pileup) are simulated for all events according to the distribution of the number of interactions observed in each bunch crossing [53]. The interactions of particles with the CMS detector are simulated using GEANT4 [54].

Object reconstruction and event selection
The data are collected using triggers based on the scalar sum of jet p T (H T ), with a requirement of H T > 900 GeV (2016) and H T > 1050 GeV (2017 and 2018). Events are reconstructed offline using a particle-flow (PF) algorithm [55] that reconstructs and identifies each individual particle (PF candidate) in an event using an optimised combination of information from the components of the CMS detector.
Jets are reconstructed by clustering the PF candidates using the anti-k T clustering algorithm [56], as implemented in the FASTJET package [57]. A distance parameter of 0.4 or 0.8 is used for standard-and large-radius jets, referred to as AK4 and AK8 jets, respectively. The jet momentum is defined as the vectorial sum of all particle momenta in the jet. To mitigate the effect of pileup, constituent charged PF candidates identified to be originating from vertices other than the primary pp interaction vertex are not used in the clustering algorithm. The primary vertex is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone, as described in Section 9.4 of Ref. [58]. For AK4 jets, an offset correction is applied to correct for remaining pileup contributions. For AK8 jets, the pileup-perparticle identification algorithm [59,60] is used to rescale the momenta of constituent neutral particles according to the probability they originated from the primary vertex. This probability is based on a local shape variable that distinguishes between collinear and soft diffuse distributions of the surrounding charged particles that are compatible with the primary vertex. For all jets, jet energy corrections are derived from simulation to bring the measured average response of jets to that of particle-level jets. In situ measurements of the momentum balance in dijet, photon+jet, Z+jet, and multijet events are used to account for any residual differences in jet energy scale and resolution between data and simulation [61,62]. Additional criteria are imposed to reject jets from spurious sources, such as electronics noise and detector malfunctions [63,64].
The identification of AK8 jets originating from two collimated b quarks (double-b tagging) is integral to the reconstruction of the H 1 . A discriminant is calculated for each jet using a double-b tagging algorithm that combines tracking and vertexing information in a multivariate approach with no strong dependence on jet mass or p T [65].
The event preselection requires two AK8 jets with p T > 170 GeV and |η| < 2.4 (so that they are within the acceptance of the tracker). If there are more than two candidate AK8 jets, the two with the largest double-b tag discriminants are selected as most likely to have originated from H 1 → bb decays. For the offline analysis, H T is defined as the scalar p T sum of all AK4 jets with p T > 40 GeV and |η| < 3.0, including AK4 jets with PF candidates clustered into AK8 jets. The H T distributions for various simulated signal and background processes are shown in Fig. 2, after implementing all preselection requirements. Since the final state contains only jets, the average signal event H T depends significantly on m SUSY , and signal events with m SUSY > 1200 GeV tend to have H T > 1500 GeV.
Additional requirements based on the expected kinematic properties of signal events are applied after the preselection. They define the kinematic event selection: 1. Both selected AK8 jets must have p T > 300 GeV and |η| < 2.4, characteristic of the jets originating from H 1 → bb decay in signal events.
2. There must be at least one AK4 jet with p T > 300 GeV and |η| < 3.0, characteristic of the quarks from squark decays in signal events. Such jets must be separated by ∆R ≡ √ (∆φ) 2 + (∆η) 2 > 1.4 from both selected AK8 jets, to avoid being constructed from the same PF candidates. 3. The event H T must exceed 1500 GeV.
Although the offline H T resolution is better than that of the trigger-level variable, the offline H T threshold is comfortably above the trigger-level H T requirements. The trigger efficiency for this analysis is measured using events collected with a single muon trigger with a muon p T threshold between 24 and 27 GeV. The efficiency for each data-taking year is nearly 100%. For the 2018 data, the |η| selection for the AK4 jets is reduced from 3.0 to 2.4 to avoid a region of the endcap electromagnetic calorimeters affected by large losses in crystal transparency, and therefore increased energy-equivalent electronics noise [66]. This change has a negligible effect on signal acceptance for all considered masses.
The fraction of signal events that satisfy the kinematic selection is essentially independent of m H 1 . It increases from about 60 to 80% as m SUSY increases from 1200 to 2000 GeV, after which it remains approximately constant.

Double-b tag based event selection
The two AK8 jets that are classified as the H 1 → bb candidates in each event are randomly assigned the labels "A" and "B". Their double-b tag discriminants define a two-dimensional (2D) parameter space, shown with simulated signal and multijet event distributions in Fig. 3. The signal events are expected to contain two H 1 → bb decays and therefore accumulate in the region where both double-b tag discriminants are large. The signal-enhanced tag region (TR) is defined as the region where the sum of the two double-b tag discriminants exceeds 1.3, illustrated by the shaded triangle in Fig. 3. Two additional regions are defined in Fig. 3 for use in the multijet background estimation and validation: the control region (CR), a multijetdominated region with negligible signal; and the validation region (VR), a more signal-like region where one of the two jets has a large double-b tag discriminant. The VR is defined sufficiently far from the TR for the signal contamination to be negligible.
About 50% of the signal events that satisfy the kinematic selection populate the TR, with vari- ation at the level of ±10% across the m H 1 and m SUSY parameter space considered. Since the multijet background is dominated by light-flavour quark and gluon initiated jets, only about 3% of these events populate the TR. For the tt, Z+jets, and W+jets backgrounds, the corresponding figures are 13, 6, and 3%, respectively.

Soft-drop mass based signal and sideband regions
In signal events, both selected AK8 jets are likely to originate from H 1 → bb decays and therefore have a jet mass close to m H 1 . The multijet background has no resonant mass peak, while the other backgrounds are only expected to exhibit peaks near the known top quark and vector bosons masses, which means that an accurate reconstruction of the jet mass is important in distinguishing signal from background. The AK8 jet masses are evaluated using the "soft-drop" algorithm [67] (with a soft-drop threshold of z cut = 0.1 and angular exponent of β = 0), in which wide-angle soft radiation is removed recursively from a jet. In signal events this algorithm achieves a relative jet mass resolution from 10% for m H 1 = 125 GeV to 20% for m H 1 = 30 GeV.
The soft-drop masses of the two AK8 jets define a 2D parameter space, shown in The event distributions for a set of signal models with different m H 1 values are shown in Fig. 5, with the signal and sideband mass regions overlaid. The peaks in the signal distributions where one or both AK8 jets have a soft-drop mass close to zero result from a selected jet originating from a single parton or one of the H 1 → bb decays lying outside the acceptance of the jet reconstruction algorithm. The latter can happen when the angular separation of the b quarks exceeds the AK8 jet distance parameter, or when the ratio of the b quark p T values is larger than 9 (such that the softer b quark would not satisfy the z cut threshold in the soft-drop algorithm). For signal models with 40 < m H 1 < 125 GeV, ≈ 50% of the events that satisfy the kinematic and TR selection fall within any of the S i . However, for m H 1 < 35 GeV the bulk of the distribution is lower in mass than S 1 , leading to a rapid decrease in signal acceptance.
The distributions of background events are also shown in Fig. 5. The majority of multijet events contain at least one AK8 jet evaluated to have a small soft-drop mass, reflecting the characteristic one-prong structure of quark and gluon jets. After applying the kinematic and TR selection criteria, approximately 5% of multijet events fall within any of the S i , with greater probability at small masses. For the vector boson and tt backgrounds the corresponding figures are 7 and 19%, respectively, concentrated in the S i corresponding to masses between the W boson and top quark masses.
For each S i there are two corresponding sideband regions, U i , used for the multijet background estimation described in Section 6. The sideband regions U 1 have a triangular form to avoid the region of very small soft-drop masses, where the density from multijet events increases sharply.

Categorisation in H T and expected yields
The selected events are classified according to three H T categories: 1500-2500, 2500-3500, and above 3500 GeV. Each H T category is divided into the 10 mass signal regions S i defined in Fig. 4, resulting in a total of 30 search regions for each data-taking year. As can be seen in Fig. 6  for TR data summed over the three data-taking years, the search region yields can be visualised through a 30-bin histogram where bins 1-10 represent the S i , in ascending order, for the first H T category. The subsequent two sets of 10 bins represent the results for the second and third H T categories. The primary background is from multijet events, estimated from data using the method described in Section 6. The expected contribution from tt events is also significant, particularly in the larger soft-drop mass regions populated by jets from hadronic top quark or W boson decays. The tt simulation is validated in a dedicated tt-enriched control region in data. In Fig. 4 this is the triangular region of the parameter space with both jet masses below 200 GeV and above the upper boundary of mass region 10. The yields from Z+jets and W+jets production are small in comparison. All expected SM backgrounds tend to exhibit small values of H T compared to signal.  Fig. 6. Although the production cross section decreases quickly with increasing m SUSY , the fraction of events in the larger-H T categories increases. Within each H T category, the distribution of events in the 10 S i bins is described by a peak with a width of about three bins, centred near the model value of m H 1 .

Multijet background estimation from data
The mass sideband regions U i form a basis for using data to estimate the multijet background. The density of the multijet background is approximately uniform within each of the 10 mass regions (spanning S i and U i for each region i illustrated in Fig. 4). Apart from U 1 , each U i is constructed to have the same area as S i such that the corresponding multijet yields, respectively denotedÛ i andŜ i , are approximately equal. The observed ratios of S i to U i yields, F i , are measured in CR data. The F i factors are found to be close to unity except for the F 1 values which are approximately 1.5.
The multijet background in the TR is estimated independently for each signal region S i : whereÛ TR i is the observed TR yield in sideband region U i after subtracting the contributions from the other simulated backgrounds. In rare cases where the predictionŜ TR i is negative, it is set equal to zero.
Since the F i factors are measured and applied in different regions of double-b tag discriminant space, any correlation between the soft-drop mass and the double-b tag discriminant of AK8 jets can bias the prediction of Eq. (1). Using a sample of data satisfying an alternative kinematic event selection with the requirement for one or more AK4 jets inverted, the variation of F i between the TR and the CR is found to be less than 10%.
The overall accuracy of the multijet estimation is assessed through closure tests. First the method is applied to simulated multijet events in the TR where, within statistical uncertainties, the predicted yields are consistent with the simulated yields for each data-taking year. Second the method is applied in the multijet-dominated VR data (defined in Fig. 3

Systematic uncertainties
The simulated events for signal and the tt, Z+jets, and W+jets backgrounds are affected by various systematic uncertainties. The efficiency for tagging (mistagging) a jet originating from two b quarks (a light-flavour quark or gluon) is corrected to match that observed in data [65]. The uncertainty in this correction corresponds to ≈ 10% in the simulated signal and background yields. The uncertainties related to the jet energy corrections are applied to the jet properties in bins of p T and η. These uncertainties affect the event H T , leading to an ≈ 4% migration of events between adjacent H T categories. The uncertainty in the soft-drop mass scale in simulation relative to data leads to a migration of events between adjacent S i and U i regions of up to 10%. The uncertainty in the simulated soft-drop mass resolution affects the widths of the simulated mass peaks. This effect is larger for signal models with small m H 1 and can reduce the S i selection efficiency by up to 20%.
The systematic uncertainties are assumed to be fully correlated among the data-taking years except for the 2016 double-b tagging uncertainties, which are assumed uncorrelated because the CMS pixel detector was upgraded prior to 2017 data-taking. Changing these correlation assumptions is found to have only a small effect on the final results. Systematic uncertainties related to integrated luminosity, pileup, PDFs, renormalisation and factorisation scales, modelling of initial-state radiation, and background cross sections were also evaluated, along with the statistical uncertainties in the simulation, and were found to make negligible contributions to the total uncertainty.
Systematic uncertainties in multijet yields arise from the systematic uncertainties in the F i factors. As described in Section 6, an uncertainty of 15% is applied to the F i in the lower two H T categories and 30% in the upper H T category, uncorrelated among different F i . Except in the lowest H T category, the total uncertainty in the multijet yield is dominated by the statistical uncertainty inÛ TR i .

Results
Binned maximum likelihood fits to the data in all 30 search regions S i for each data-taking year are carried out under background-only and signal+background hypotheses. The corresponding sideband regions U i are fitted simultaneously, thereby constraining the multijet contributions to the search region yields through Eq. (1). The likelihood functions are defined through the product of 90 × 2 Poisson distributions [68], one for each search region and one for each sideband region, with additional constraint terms for the "nuisance" parameters that account for the systematic uncertainties summarised in Section 7. Figure 8 compares the result of the background-only fit to the yields in the search regions for the combination of 2016, 2017, and 2018 data. There is no evidence for deviations of the data from the fitted background. The values and uncertainties of most nuisance parameters are unchanged in the fit, but the ones corresponding to the F i are constrained through Eq. (1) when the yieldsŜ TR i andÛ TR i are sufficiently large.
Signal+background fits are used to set 95% confidence level (CL) upper limits on the product σB 2 for the mass points in the benchmark signal model. The limits are set using the modified frequentist CL s criterion [69,70], with the profile likelihood ratio as test statistic [68]. The observed and expected 95% CL upper limits on σB 2 are shown in Fig. 9, as functions of m H 1 for constant m SUSY . The upper limits are weaker for models with m H 1 < 35 GeV, for which the signal-event distribution in the 2D soft-drop mass plane peaks outside the signal regions. The limits have no significant dependence on m SUSY for models with m SUSY > 2000 GeV, whose signal events mostly populate the upper H T category (as shown in Fig. 6).
The σB 2 upper limits are used in conjunction with the theoretical σ and B values from Section 2 to exclude ranges of masses in m H 1 and m SUSY in the benchmark model. The observed 95% CL upper limits on r, the ratio of measured and theoretical values of σB 2 , are shown in Fig. 10, with the corresponding exclusion contours at r = 1. Masses 1200 < m SUSY < 2500 GeV are excluded within the range 40 < m H 1 < 120 GeV. Expected exclusion contours for the backgroundonly scenario agree within one standard deviation with the observed contours. In the region 110 < m H 1 < 125 GeV, B starts to decrease more quickly (as shown in Table 2), leading to a corresponding reduction in sensitivity. Most of the sensitivity at large m SUSY comes from the H T > 3500 GeV region, where the statistical uncertainties in the observed yields are dominant over systematic uncertainties. This search does not explore the region outside of that shown in Fig. 10.
To aid reinterpretation of the search by reducing the model-dependence, limits evaluated using only the upper H T category are presented in Appendix A. Tabulated results are provided in the HEPData record for this analysis [71].     [17] CMS Collaboration, "Search for higgsino pair production in pp collisions at √ s = 13 TeV in final states with large missing transverse momentum and two Higgs bosons decaying via H → bb", Phys. Rev. D 97 (2018) 032007, doi:10.1103/PhysRevD.97.032007, arXiv:1709.04896.

A Simplified analysis for reinterpretation
To aid reinterpretation of the search, a simplified analysis is performed using only the 10 search regions in the upper H T category. The value A kin is defined as the product of acceptance and efficiency for a signal event to satisfy the kinematic selection (defined in Section 5) and the H T > 3500 GeV requirement. The value of A kin is common among all 10 search regions in the simplified analysis, and is quoted for the benchmark signal model in Table A.1. Upper limits on the product σ B 2 A kin as a function of m H 1 are set in Fig. A.1, from which σB 2 limits for different signal models can be derived through division by the appropriate value of A kin . Since the upper H T category provides most of the sensitivity for m SUSY > 2000 GeV in the nominal analysis, the σB 2 upper limits in this region are not much weaker in the simplified analysis. This is not the case in the region m SUSY < 2000 GeV, where the lower H T categories become important.  The observed and expected 95% CL upper limit on the product of σB 2 and A kin , the kinematic acceptance and efficiency for the H T > 3500 GeV region, as a function of m H 1 . The results are independent of m SUSY within 10% in the range 1600 < m SUSY < 2800 GeV. The solid and dashed black lines indicate the observed and median expected limits, respectively. The inner (green) band and the outer (yellow) band indicate the regions containing 68 and 95%, respectively, of the distribution of limits expected under the background-only hypothesis.
The double-b tag and mass region selections are not considered in A kin . This is done for simplicity, and because the fraction of events satisfying these selections is not found to be strongly model-dependent (except for the dependence on m H 1 , which is accounted for explicitly in Fig. A.1). For the benchmark model, this fraction is found to be independent of m SUSY within 10% in the region 1600 < m SUSY < 2800 GeV and 35 < m H 1 < 125 GeV. This approximate independence does not hold for models with m SUSY < 1600 GeV, where the H 1 p T distribution has substantial contributions below the p T necessary for the H 1 → bb decay products to be merged in a single AK8 jet. Only models with typical bb angular separation ∆R < 0.8 should be considered for reinterpretation.