Search for natural and split supersymmetry in proton-proton collisions at $\sqrt{s} =$ 13 TeV in final states with jets and missing transverse momentum

A search for supersymmetry (SUSY) is performed in final states comprising one or more jets and missing transverse momentum using data from proton-proton collisions at a centre-of-mass energy of 13 TeV. The data were recorded with the CMS detector at the CERN LHC in 2016 and correspond to an integrated luminosity of 35.9 fb$^{-1}$. The number of signal events is found to agree with the expected background yields from standard model processes. The results are interpreted in the context of simplified models of SUSY that assume the production of gluino or squark pairs and their prompt decay to quarks and the lightest neutralino. The masses of bottom, top, and mass-degenerate light-flavour squarks are probed up to 1050, 1000, and 1325 GeV, respectively. The gluino mass is probed up to 1900, 1650, and 1650 GeV when the gluino decays via virtual states of the aforementioned squarks. The strongest mass bounds on the neutralinos from gluino and squark decays are 1150 and 575 GeV, respectively. The search also provides sensitivity to simplified models inspired by split SUSY that involve the production and decay of long-lived gluinos. Values of the proper decay length $c\tau$ from $10^{-3}$ to $10^{5}$ mm are considered, as well as a metastable gluino scenario. Gluino masses up to 1750 and 900 GeV are probed for $c\tau =$ 1 mm and for the metastable state, respectively. The sensitivity is moderately dependent on model assumptions for $c\tau \gtrsim 1$ m. The search provides coverage of the $c\tau$ parameter space for models involving long-lived gluinos that is complementary to existing techniques at the LHC.


Introduction
Supersymmetry (SUSY) [1][2][3][4] is an extension of the standard model (SM) of particle physics that introduces at least one bosonic (fermionic) superpartner for each fermionic (bosonic) SM particle, where the superpartner differs in its spin from its SM counterpart by one half unit. Supersymmetry offers a potential solution to the hierarchy problem [5,6], predicts unification of the gauge couplings at high energy [7][8][9], and provides a candidate for dark matter (DM). Under the assumption of R-parity [10] conservation, SUSY particles are expected to be produced in pairs at the CERN LHC and to decay to the stable, lightest SUSY particle (LSP). The LSP is assumed to be the neutralino χ 0 1 , a weakly interacting massive particle and a viable DM candidate [11,12]. So-called natural SUSY models, which invoke only a minimal fine tuning of the bare Higgs boson mass parameter, require only the gluino, third-generation squarks, and a higgsino-like χ 0 1 to have masses at or near the electroweak (EW) scale [13]. The interest in natural models is motivated by the discovery of a low-mass Higgs boson [14][15][16][17][18][19]. The characteristic signature of natural SUSY production at the LHC is a final state containing an abundance of jets originating from the hadronization of heavy-flavour quarks and significant missing transverse momentum p miss T . Split supersymmetry [20,21] does not address the hierarchy problem-in contrast to natural SUSY models-but preserves the appealing aspects of gauge coupling unification and a DM candidate. In such a model, only the fermionic superpartners, and a finely tuned scalar Higgs boson, may be realized at a mass scale that is kinematically accessible at the LHC. All other SUSY particles are assumed to be ultraheavy. Hence, within split SUSY models, the gluino decay is suppressed because of the highly virtual squark states. For gluino lifetimes beyond a picosecond, the gluino hadronizes and forms a bound colour-singlet state containing the gluino and quarks or gluons [22], known as an R-hadron, before eventually decaying to a quarkantiquark pair and the χ 0 1 . The long-lived gluino can lead to final states with significant p miss T from the undetected χ 0 1 particles and to jets with vertices located a significant distance (i.e. displaced) from the luminous region of the proton beams. A metastable gluino, with a decay length significantly beyond the scale of the CMS detector, can escape undetected. This paper presents a search for new-physics processes in final states with one or more energetic jets and significant p miss T . The search is performed with a sample of proton-proton (pp) collision data at a centre-of-mass energy of 13 TeV recorded by the CMS experiment in 2016. The analysed data sample corresponds to an integrated luminosity of 35.9 ± 0.9 fb −1 [23]. Earlier searches using the same technique have been performed in pp collisions at √ s = 7, 8, and 13 TeV by the CMS Collaboration [24][25][26][27][28][29]. The data set analysed in this analysis is a factor of 16 larger than that presented in Ref. [29]. The search strategy aims to provide sensitivity to a broad range of SUSY-inspired models that predict the existence of a DM candidate, and the search is used to constrain the parameter spaces of a number of simplified SUSY models [30][31][32]. The overwhelmingly dominant background for a new-physics search in all-jet final states resulting from pp collisions is the production of multijet events via the strong interaction, a manifestation of quantum chromodynamics (QCD). Several dedicated variables are employed to suppress the multijet background to a negligible level while maintaining low kinematical thresholds and high experimental acceptance for final states characterized by the presence of significant p miss T . Signal extraction is performed using additional kinematical variables, namely the number of jets, the number of jets identified as originating from bottom quarks, and the scalar and vector sums of the jet transverse momenta. The ATLAS and CMS Collaborations have performed similar searches in all-jet final states at √ s = 13 TeV, of which those providing the tightest constraints are described in Refs. [33][34][35]. This search does not employ specialized reconstruction techniques [36][37][38][39][40][41][42][43][44][45][46] that target long-lived gluinos. This paper is organized as follows. Section 2 describes the CMS apparatus and the event reconstruction algorithms. Section 3 summarizes the selection criteria used to identify and categorize signal events and samples of control data. Section 4 outlines the various software packages used to generate the samples of simulated events. Sections 5 and 6 describe the methods used to estimate the background contributions from SM processes. The results and interpretations are described in Sections 7 and 8, respectively, and summarized in Section 9.

The CMS detector and event reconstruction
The central feature of the CMS detector is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematical variables, can be found in Ref. [47].
Events of interest are selected using a two-tiered trigger system [48]. 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 time interval of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to less than 1 kHz before data storage. The trigger logic used by this search is summarized in Section 3.
The particle-flow (PF) event algorithm [49] reconstructs and identifies each individual particle with an optimized combination of information from the various elements of the CMS detector. In this process, the identification of the particle type (photon, electron, muon, charged hadron, neutral hadron) plays an important role in the determination of the particle direction and energy. The energy of photons [50] is directly obtained from the ECAL measurement, corrected for zero-suppression effects. The energy of electrons [51] is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL measurement, and the energy sum of all bremsstrahlung photons spatially compatible with originating from the electron track. The energy of muons [52] is obtained from the curvature of the corresponding track. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for zero-suppression effects and for the response function of the calorimeters to hadronic showers. Finally, the energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energy. The reconstruction techniques used by this search are not specialized to target specific experimental signatures (such as displaced jets). The physics objects used in this search are defined below and are summarized in Table 1. In the case of photons and leptons, further details can be found in Ref. [29] and references therein.
The reconstructed vertex with the largest value of summed physics object p 2 T is taken to be the primary pp interaction vertex (PV). The physics objects considered are those returned by a jet finding algorithm [53,54] applied to all charged particle tracks associated with the vertex, and the associated p miss T , taken as the negative vector sum of the p T of those physics objects. Charged particle tracks associated with vertices from additional pp interactions within the same or nearby bunch crossings (pileup) are not considered by the PF algorithm as part Table 1: Summary of the physics object acceptances, the baseline event selection, the signal and control regions, and the event categorization schemas. The nominal categorization schema is defined in full in Appendix A.

Physics object acceptances
Sidebands to signal region: H miss T /p miss T > 1.25 and/or ∆φ * min < 0.5 of the global event reconstruction. The energy deposit associated with each physics object is corrected to account for contributions from neutral particles originating from pileup interactions [55].
Samples of signal events and control data are defined, respectively, by the absence or presence of photons and leptons that are isolated from other activity in the event. Photons are required to be isolated [50] within a cone around the photon trajectory defined by the radius ∆R = √ (∆φ) 2 + (∆η) 2 = 0.3, where ∆φ and ∆η represent differences in the azimuthal angle (radians) and pseudorapidity. Isolation for an electron or muon is a relative quantity, I rel , defined as the scalar p T sum of all PF particle candidates within a cone around its trajectory, divided by the lepton p T . The cone radius is dependent on the lepton p T , ∆R = min[max(0.05, 10 GeV/p T ), 0.2], to maintain high efficiency for semileptonic decays of Lorentz-boosted top quarks [56]. Isolated electrons and muons are required to satisfy I rel < 0.1 and 0.2, respectively. Electron and muon candidates that fail any of the aforementioned requirements, as well as charged-hadron candidates from hadronically decaying tau leptons, are collectively labelled as single isolated tracks (SITs) if the scalar p T sum of additional tracks associated with the PV within a cone ∆R < 0.3 around the track trajectory, relative to the track p T , satisfies I track < 0.1. All isolation variables exclude the contributions from the physics object itself and pileup interactions [50][51][52]. The experimental acceptances for photons, electrons, muons, and SITs are defined by the transverse momentum requirements p T > 25, 10, 10, and 10 GeV, respectively, and the pseudorapidity requirement |η| < 2.5.
Jets are reconstructed from the PF particle candidates, clustered by the anti-k T algorithm [53,54] with a distance parameter of 0.4. In this process, the raw jet energy is obtained from the sum of the particle candidate energies, and the raw jet momentum by the vectorial sum of the particle candidate momenta, which results in a nonzero jet mass. An offset correction is applied to jet energies to take into account the contributions from neutral particles produced in pileup interactions [55,57]. The raw jet energies are then corrected to establish a relative uniform response of the calorimeter in η and a calibrated absolute response in p T . Jet energy corrections are derived from simulation, and are confirmed with in situ measurements of the energy balance in dijet, multijet, γ+jets, and leptonically decaying Z+jets events [58]. Jets are required to satisfy p T > 40 GeV and |η| < 2.4. Jets are also subjected to a standard set of identification criteria [59] that require each jet to contain at least two particle candidates and at least one charged particle track, and the energy fraction f h ± attributed to charged-hadron particle candidates is required to be nonzero.
Jets can be identified as originating from b quarks using the combined secondary vertex (CSVv2) algorithm [60]. Data samples are used to measure the b tagging efficiency, which is the probability to correctly identify jets originating from b quarks, as well as the mistag probabilities for jets that originate from light-flavour (LF) partons (u, d, s quarks or gluon) or a charm quark. A working point is employed that yields a b tagging efficiency of ≈69% for jets with p T > 30 GeV from tt events, and charm and LF mistag probabilities of ≈18 and ≈1%, respectively, for multijet events.
Finally, the most accurate estimator for p miss T is defined as the projection on the plane perpendicular to the beams of the negative vector sum of the momenta of all PF particle candidates in an event. Its magnitude is referred to as p miss T .

Event selection and categorization
A baseline set of event selection criteria, described in Section 3.1, is used as a basis for all data samples used in this search. Two additional requirements, described in Section 3.2, are employed to define a sample of signal events, labelled henceforth as the signal region (SR). The categorization of signal events and the background composition are described in Sections 3.3 and 3.4, respectively. Three independent control regions (CRs), comprising large samples of event data, are defined by the selection criteria described in Section 3.5. All selection criteria are summarized in Table 1.

Baseline selections
Events containing isolated photons, electrons and muons, or SITs that satisfy the requirements summarized in Table 1 are vetoed. The aforementioned vetoes are employed to select all-jet final states, suppress SM processes that produce final states containing neutrinos, and reduce backgrounds from misreconstructed or nonisolated leptons as well as single-prong hadronic decays of τ leptons.
Beam halo, spurious jet-like features originating from isolated noise patterns in the calorimeter systems, detector inefficiencies, and reconstruction failures can all lead to large values of p miss T . Such events are rejected with high efficiency using dedicated vetoes [61,62]. Events are vetoed if any jet fails the identification criteria described in Section 2. Further, f h ± for the highest p T jet of the event, j 1 , is required to satisfy 0.1 < f j 1 h ± < 0.95 to further suppress beam halo and rare reconstruction failures.
The highest p T jet in the event is required to satisfy p variable. An additional veto is employed to deal with the circumstance in which several jets with transverse momentum below the p T thresholds and collinear in φ can result in significant H miss T relative to p miss T , the latter of which is less sensitive to jet thresholds. This type of event topology, which is typical of multijet events, is suppressed while maintaining high efficiency for new-physics processes with significant p miss T by requiring H miss T /p miss T < 1.25.

Signal region
The multijet background dominates over all other SM backgrounds following the application of the baseline event selection criteria. The multijet background is suppressed to a negligible level through the application of two dedicated variables that provide strong discrimination between multijet events with p miss T resulting from instrumental sources, such as jet energy mismeasurements, and new-physics processes that involve the production of weakly interacting particles that escape detection.
The first variable, α T [24, 63], is designed to be intrinsically robust against jet energy mismeasurements. In its simplest form, the α T variable is defined as T (1 − cos φ j 1 ,j 2 ) and φ j 1 ,j 2 is defined as the azimuthal angle between jets j 1 and j 2 . In the absence of jet energy mismeasurements, and in the limit for which the E T of each jet is large compared with its mass, a well-measured dijet event with E j 1 T = E j 2 T and back-to-back jets (φ j 1 ,j 2 = π) yields an α T value of 0.5. In the presence of a jet energy mismeasurement, E j 1 T > E j 2 T and α T < 0.5. Values significantly greater than 0.5 can be observed when the two jets are not back-to-back and recoil against p miss T from weakly interacting particles that escape the detector. The definition of the α T variable can be generalized for events with two or more jets, as described in Ref. [24]. Multijet events populate the region α T 0.5 and the α T distribution is characterized by a sharp edge at 0.5, beyond which the multijet event yield falls by several orders of magnitude. The SM backgrounds that involve the production of neutrinos result in a long tail in α T beyond values of 0.5. A H T -dependent α T threshold that decreases from 0.65 at low H T to 0.52 at high H T within the range 200 < H T < 900 GeV is employed to maintain an approximately constant rejection power against the multijet background.
The second variable, known as ∆φ * min , considers the minimum azimuthal angular separation between each jet in the event and the vector p T sum of all other jets in the event. Multijet events typically populate the region ∆φ * min ≈ 0 while events with genuine p miss T can have values up to ∆φ * min = π. The requirement ∆φ * min > 0.5 is sufficient to reduce significantly the multijet background, including rare contributions from energetic multijet events that yield both high jet multiplicities and significant p miss T because of high-multiplicity neutrino production in semileptonic heavy-flavour decays. For events that satisfy n jet = 1, a small modification to the ∆φ * min variable is utilized that considers any additional jets with 25 < p T < 40 GeV that are outside the nominal experimental acceptance (∆φ * 25 min > 0.5). The requirements on α T and ∆φ * min , summarized in Table 1, suppress the expected contribution from multijet events to the sub-percent level with respect to the total expected background counts from all other SM processes. For the region H T > 900 GeV, the necessary control of the multijet background is achieved solely with the ∆φ * min and ∆φ * 25 min variables. These requirements complete the definition of the SR.
Signal events are recorded with a number of trigger algorithms. Events with n jet ≥ 2 must satisfy thresholds on both H T and α T that are looser than those used to define the SR. Highactivity events that satisfy H T > 900 GeV are also recorded. Finally, a trigger condition that requires H miss T > 120 GeV, p miss T > 120 GeV, and a single jet with p T > 20 GeV and |η| < 5.2 is also used to efficiently record signal events for all categories of the SR, including those that satisfy n jet ≥ 1. The combined performance of these trigger algorithms yields high efficiencies, as determined from samples of CR data enriched in vector boson + jets and tt events. The efficiencies are primarily H T -dependent and range from 97.4-97.9% (200 < H T < 600 GeV) to 100% (H T > 600 GeV) with statistical and systematic uncertainties at the percent level. Trigger efficiencies for a range of benchmark signal models are typically comparable or higher (≈100%).

Event categorization
Signal events are categorized into 27 discrete topologies according to n jet and the number of b-tagged jets n b . Events are further binned according to the energy sums H T and H miss T . The binning schema is determined primarily by the statistical power of the µ+jets and µµ+jets CRs.
Seven bins in n jet are considered, as summarized in Table 1. Events that contain only a single jet within the experimental acceptance (n jet = 1) are labelled as "monojet". Events containing two or more jets are categorized according to the second-highest jet p T . Events that satisfy n jet ≥ 2 with only the highest p T jet satisfying p T > 100 GeV are labelled as "asymmetric". Events for which the second-highest jet p T also satisfies p T > 100 GeV are labelled as "symmetric" and are categorized according to n jet (2, 3, 4, 5, and ≥6). The symmetric topology targets the pair production of SUSY particles and their prompt cascade decays, while the monojet and asymmetric topologies preferentially target models with a compressed mass spectrum and long-lived SUSY particles.
Events are also categorized according to n b (0, 1, 2, 3, ≥4), where n b is bounded from above by n jet and the choice of categorization is dependent on n jet . Higher n b multiplicities target the production of third-generation squarks.
The nominal binning schema for H T is defined as follows: four bounded bins that satisfy 200-400, 400-600, 600-900, and 900-1200 GeV, and a final open bin H T > 1200 GeV. This schema is adapted per (n jet , n b ) category as follows: only the region H T > 400 GeV is considered for events that satisfy n jet ≥ 4, and bins at high H T are merged with lower-H T bins to satisfy a threshold of at least four events in the corresponding bins of the CRs.
The H miss T variable is used to further categorize events according to three bounded bins that In total, there are 254 bins in the SR. An alternate, simplified binning schema is also provided in which events are categorized according to eight topologies defined in terms of n jet and n b . For each topology, event yields are integrated over the full available H T range and categorized according to the four nominal H miss T bins defined above. This schema has 32 bins that are exclusive, contiguous, and provide a complete coverage of the SR. The SM background estimates are obtained from the same likelihood model as the one used to determine the nominal result.

Background composition
Following the application of the SR selection criteria, the multijet background is reduced to a negligible level. The remaining background contributions are dominated by processes that involve the production of high-p T neutrinos in the final state. The associated production of jets and a Z boson that decays to νν dominates the background contributions for events containing low numbers of jets and b-tagged jets. The Z(→ νν)+jets background is irreducible. The associated production of jets and a W boson that decays to ν ( = e, µ, τ) is also a significant background in the same phase space. The production and semileptonic decay of top quark-antiquark pairs (tt) becomes the dominant background process for events containing high numbers of jets or b-tagged jets. Events that contain the leptonic decay of a W boson are typically rejected by the vetoes that identify the presence of leptons or single isolated tracks. If the lepton is outside the experimental acceptance, or is not identified or isolated, then the event is not vetoed and the aforementioned processes lead to what is collectively known as the "lost lepton" ( lost ) background. Residual contributions from other SM processes are also considered, such as single top quark production; WW, WZ, ZZ (diboson) production; and the associated production of tt and a boson (ttW, ttZ, ttγ, and ttH).

Control regions
Topological and kinematical requirements, summarized in Table 1, ensure that the samples of CR data are enriched in the same or similar SM processes that populate the SR, as well as being depleted in contributions from SUSY processes (signal contamination).
Three sidebands to the SR comprising multijet-enriched event samples are defined by , and both 1.25 < H miss T /p miss T < 3.0 and 0.2 < ∆φ * min < 0.5 (C). Events are categorized according to n jet and H T , identically to the SR. Events are recorded with the signal triggers described above.
Two additional CRs comprising µ+jets and µµ+jets event samples are defined by the application of the baseline selections and requirements on isolated, central, high-p T muons. Tighter isolation requirements for the muons are applied with respect to those indicated in Table 1. A trigger condition that requires an isolated muon with p T > 24 GeV and |η| < 2.1 is used to record the µ+jets and µµ+jets event samples with efficiencies of ≈90 and ≈99%, respectively. For both samples, no requirements on α T or ∆φ * min are imposed. The kinematical properties of events in the µ+jets and µµ+jets CRs and SR are comparable once the muon or dimuon system is ignored in the calculation of event-level quantities such as H T and H miss T . Events in both samples are categorized according to n jet , H T , and n b , with counts integrated over H miss T . The n jet categorization is identical to the SR. Background predictions are determined using up to eleven bins in H T that are then aggregated to match the H T binning schema used by the SR. The n b categorization for µ+jets events is identical to the SR, whereas µµ+jets events are subdivided according to n b = 0 and n b ≥ 1. Differences in the binning schemas between the SR and CRs are accounted for in the background estimation methods through simulation-based templates, the modelling of which is validated against control data.
The µ+jets event sample is enriched in events from W(→ µν)+jets and tt production, as well as other SM processes (e.g. single top quark and diboson production), that are manifest in the SR as the lost backgrounds. Each event is required to contain a single isolated muon with p T > 30 GeV and |η| < 2.1 to satisfy trigger conditions, and is well separated from each jet j i in the event according to ∆R(µ, The µµ+jets sample is enriched in Z → µ + µ − events that have similar acceptance and kinematical properties to Z(→ νν)+jets events when the muons are ignored. The sample uses selection criteria similar to the µ+jets sample, but requires two oppositely charged, isolated muons that both satisfy p T > 30 GeV, |η| < 2.1, and ∆R(µ 1,2 , j i ) > 0.5. The muons are also required to have a dilepton invariant mass m µµ within a ±25 GeV window around the mass of the Z boson [12].

Monte Carlo simulation
The search relies on several samples of simulated events, produced with Monte Carlo (MC) generator programs, to aid the estimation of SM backgrounds and evaluate potential signal contributions.
The MADGRAPH5 aMC@NLO 2.2.2 [64] event generator is used at leading-order (LO) accuracy to produce samples of W+jets, Z+jets, tt, and multijet events. Up to three or four additional partons are included in the matrix-element calculation for tt and vector boson production, respectively. Simulated W+jets and Z+jets events are weighted according to the true vector boson p T to account for the effect of missing next-to-leading-order (NLO) QCD and EW terms in the matrix-element calculation [64,65], according to the method described in Ref. [66]. Within the range of vector boson p T that can be probed by this search, the QCD and EW corrections [65] are largest, ≈40% and ≈15%, at low and high values of boson p T , respectively. Simulated tt events are weighted to improve the description of jets arising from initial-state radiation (ISR) [67]. The weights vary from 0.92 to 0.51 depending on the number of jets (1-6) from ISR, with an uncertainty of one half the deviation from unity. The MADGRAPH5 aMC@NLO generator is used at NLO accuracy to generate samples of s-channel production of single top quark, as well as ttW and ttZ events. The NLO POWHEG v2 [68,69] generator is used to describe the tand Wtchannel production of events containing single top quarks, as well as ttH events. The PYTHIA 8.205 [70] program is used to generate diboson (WW, WZ, ZZ) production.
Event samples for signal models involving the production of gluino or squark pairs, in association with up to two additional partons, are generated at LO with MADGRAPH5 aMC@NLO, and the decay of the SUSY particles is performed with the PYTHIA program. The NNPDF3.0 LO and NNPDF3.0 NLO [71] parton distribution functions (PDFs) are used, respectively, with the LO and NLO generators described above.
The simulated samples for SM processes are normalized according to production cross sections that are calculated with NLO and next-to-NLO precision [64,69,[72][73][74][75][76]. The production cross sections for pairs of gluinos or squarks are determined at NLO plus next-to-leading-logarithm (NLL) precision [77][78][79][80][81][82]. All other SUSY particles, apart from the χ 0 1 , are assumed to be heavy and decoupled from the interaction. Uncertainties in the cross sections are determined from different choices of PDF sets, and factorization and renormalization scales (µ F and µ R ), according to the prescription in Ref. [82]. The PYTHIA program with the CUETP8M1 tune [83,84] is used to describe parton showering and hadronization for all simulated samples.
The RHADRONS package within the PYTHIA 8.205 program is used to describe the formation of R-hadrons through the hadronization of gluinos [22, 85,86]. The hadronization process, steered according to the default parameter settings of the RHADRONS package, predominantly yields meson-like ( gqq) and baryon-like ( gqqq) states, as well as glueball-like ( gg) states with a probability P gg = 10%, where g, g, q, and q represent a gluino, gluon, quark, and antiquark, respectively. The gluino is assumed to undergo a three-body decay, to a qq pair and the χ 0 1 , according to its proper decay length cτ 0 that is a parameter of the simplified model [87]. Studies with alternative values for parameters that influence the hadronization of the gluino, such as P gg = 50%, indicate a minimal influence on the event topology and kinematical variables for the models considered in this paper. Further, the model-dependent interactions of R-hadrons with the detector material are not considered by default, as studies demonstrate that the sensitivity of this search is only moderately dependent on these interactions, as discussed in Section 8.
The description of the detector response is implemented using the GEANT4 [88] package for all simulated SM processes. Scale factors are applied to simulated event samples that correct for differences with respect to data in the b tagging efficiency and mistag probabilities. The scale factors have typical values of ≈0.95-1.00 and ≈1.00-1.20, respectively, for a jet p T range of 40-600 GeV [60]. All remaining signal models rely on the CMS fast simulation package [89] that provides a description that is consistent with GEANT4 following the application of near-unity corrections for the differences in the b tagging efficiency and mistag probabilities, as well as corrections for the differences in the modelling of the H miss T distribution. To model the effects of pileup, all simulated events are generated with a nominal distribution of pp interactions per bunch crossing and then weighted to match the pileup distribution as measured in data.

Nonmultijet background estimation
The lost and Z(→ νν)+jets backgrounds, collectively labelled henceforth as the nonmultijet backgrounds, are estimated from data samples in CRs and transfer factors R determined from the ratios of expected counts obtained from simulation: where R lost and R Z→νν are the transfer factors that act as multiplier terms on the event counts N µ+jets data and N µµ+jets data observed in each (n jet , H T , n b ) bin of, respectively, the µ+jets and µµ+jets CRs to estimate the lost or Z(→ νν)+jets background counts N lost pred and N Z→νν pred in the corresponding (n jet , H T , n b , H miss T ) bins of the SR. Several sources of uncertainty in the transfer factors are evaluated. In addition to statistical uncertainties arising from finite-size simulated event samples, the most relevant systematic effects are discussed below.  The uncertainties from known theoretical and experimental sources are propagated through to the transfer factors to ascertain the magnitude of variations related to the following: the jet energy scale, the scale factors related to the b tagging efficiency and mistag probabilities, the efficiency to trigger on and identify, or veto, well-reconstructed isolated leptons, the PDFs [90], µ F and µ R , and the modelling of jets from ISR produced in association with tt [67]. Uncertainties of 100% in both the NLO QCD and EW corrections to the W+jets and Z+jets simulated samples are also considered. A 5% uncertainty in the total inelastic cross section [91] is assumed and propagated through to the weighting procedure to account for differences between the data and simulation in the pileup distributions. Uncertainties in the signal trigger efficiency measurements are also propagated to the transfer factors. The effects of the aforementioned systematic uncertainties are summarized in Table 2, in terms of representative ranges. Each source of uncertainty is assumed to vary with a fully correlated behaviour across the full phase space of the SR and CRs.
Sources of additional uncertainties are determined from closure tests performed using control data that aim to identify n jet -or H T -dependent sources of systematic bias arising from extrapolations in kinematical variables covered by the transfer factors. Several sets of tests are performed. The accuracy of the modelling of the efficiencies of both the α T and ∆φ * min requirements is estimated from both the µ+jets and µµ+jets samples. The effects of W boson polarization are probed by using µ+jets events with a positively charged muon to predict those containing a negatively charged muon. Finally, the efficiency of the single isolated track veto is also probed using a sample of µ+jets events. The uncertainties are summarized in Table 2.
The simulation modelling of the n b distributions for the Z(→ νν)+jets background in the region n b ≥ 1 is evaluated through a binned maximum-likelihood fit to the observed n b distributions in data in each (n jet , H T ) bin of the µµ+jets CR. Additional checks are performed in µµ+jets samples that are enriched in mistagged jets that originate from LF partons or charm quarks, or the genuine tags of b quarks from gluon splitting, through the use of loose and tight working points of the b tagging algorithm, respectively. No tests reveal evidence of significant bias in the simulation modelling of the n b distribution.
Finally, the modelling of the H miss T distribution in simulated events is compared to the distributions observed in µ+jets and µµ+jets control data, and inspected for trends, by assuming a linear behaviour of the ratio of observed and simulated counts as a function of H miss T . Linear fits are performed independently for each n jet category while integrating event counts over n b and H T , and then repeated for each H T bin while integrating event counts over n jet and n b . Systematic uncertainties are determined from any nonclosure between data and simulation as a function of n jet and are assumed to be correlated in H T (and n b ), and vice versa. The uncertainties can be as large as ≈50% in the most sensitive H miss T bins.

Multijet background estimation
The multijet background is estimated using the three data sidebands defined in Section 3.5. Events in each sideband are categorized according to n jet and H T . The event counts in data are corrected to account for contamination from nonmultijet SM processes, such as vector boson and tt production, as well as the residual contributions from other SM processes. The nonmultijet processes are estimated from the µ+jets and µµ+jets CRs, following a procedure similar to the one described in Section 5. The corrected counts are assumed to arise solely from multijet production. For each sideband, a transfer factor per (n jet , H T ) bin is obtained from simulation, defined as the ratio of the number of multijet events that satisfies the SR requirements to the number that satisfies the sideband requirement. Estimates of the multijet background per (n jet , H T ) bin are obtained per sideband from the product of the transfer factors and the corrected data counts.
The final estimate per (n jet , H T ) bin is a weighted sum of the three estimates. The multijet background is found to be small, typically at the percent level, relative to the sum of all nonmultijet backgrounds in all (n jet , n b ) bins of the SR. The H miss T /p miss T and ∆φ * min variables that are used to define the sidebands are determined to be only weakly correlated for multijet events, and the estimates from each sideband are assumed to be uncorrelated. Statistical uncertainties associated with the finite event counts in data and simulated event samples, as large as ≈100%, are propagated to each estimate. Uncertainties as large as ≈20% in the estimates of nonmultijet contamination are also propagated to the corrected events. Any differences between the three estimates per (n jet , H T ) bin are adequately covered by systematic uncertainties of 100%, which are assumed to be uncorrelated across (n jet , H T ) bins.
A model is assumed to determine the estimates as a function of n b and H miss T . The distribution of multijet events as a function of n b and H miss T per (n jet , H T ) bin is assumed to be identical to the distribution expected for the nonmultijet backgrounds. This assumption is based on simulation-based studies and is a valid simplification given the magnitude of the multijet background relative to the sum of all other SM backgrounds, as well as the magnitude of the statistical and systematic uncertainties in the estimates described above.

Results
A likelihood model is used to obtain the SM expectations in the SR and each CR, as well as to test for the presence of new-physics signals. The observed event count in each bin, defined in terms of the n jet , n b , H T , and H miss T variables, is modelled as a Poisson-distributed variable around the SM expectation and a potential signal contribution (assumed to be zero in the following discussion). The expected event counts from nonmultijet processes in the SR are related to those in the µ+jets and µµ+jets CRs via simulation-based transfer factors, as described in Section 5. The systematic uncertainties in the nonmultijet estimates, summarized in Table 2, are accommodated in the likelihood model as nuisance parameters, the measurements of which are assumed to follow a log-normal distribution. In the case of the modelling of the H miss T distribution, alternative templates are used to describe the uncertainties in the modelling and a vertical template morphing schema [29,92] is used to interpolate between the nominal and alternative templates. The multijet background estimates, determined using the method described in Section 6, are also included in the likelihood model. Figures 1 and 2 summarize the binned counts of signal events and the corresponding SM expectations as determined from a "CR-only" fit that uses only the data counts in the µ+jets and µµ+jets control regions to constrain the model parameters related to the nonmultijet backgrounds. The uncertainties in the SM expectations reflect both statistical and systematic components. The multijet background estimates are determined independently and included in the SM expectations. The fit does not consider the event counts in the signal region.
Hypothesis testing with regards to a potential signal contribution is performed by considering a full fit to the event counts in the SR and CRs. No significant deviation is observed between the predictions and data in the SR and CRs, and the data counts appear to be adequately modelled by the SM expectations with no significant kinematical patterns. Event counts, SM background estimates, and the associated correlation matrix are also determined using the simplified 32-bin schema, which can be found in Appendix A.

Interpretations
The search result is used to constrain the parameter spaces of simplified SUSY models [30-32]. Interpretations are provided for nine unique model families, as summarized in Table 3. Each family of models realizes a unique production and decay mode. The model parameters are the masses of the parent gluino (m g ) or bottom, top, and LF (m b 1 , m t 1 , m q ) squark, also collectively labelled as m SUSY , and the χ 0 1 (m χ 0 1 ). Two scenarios are considered for LF squarks: one with an eightfold mass degeneracy for q L and q R with q = { u, d, s, c} and the other with just a single light squark ( u L ). All other SUSY particles are assumed to be too heavy to be produced directly.
Gluinos are assumed to undergo prompt three-body decays via highly virtual squarks. In the case of split SUSY models (T1qqqqLL), the gluino is assumed to be long-lived with proper decay lengths in the range 10 −3 < cτ 0 < 10 5 mm. A scenario involving a metastable gluino with cτ 0 = 10 18 mm is also considered.
Under the signal+background hypothesis, and in the presence of a nonzero signal contribution, a modified frequentist approach is used to determine observed upper limits (ULs) at 95% confidence level (CL) on the cross section σ UL to produce pairs of SUSY particles as a function of

Model family Production and decay Additional assumptions
Production and prompt decay of squark pairs

Production and prompt decay of gluino pairs
pp → g g, g → q q * → qq χ 0 1 m q m g Production and decay of long-lived gluino pairs T1qqqqLL pp → g g, g → q q * → qq χ 0 1 m q m g , 10 −3 < cτ 0 < 10 5 mm or metastable m SUSY , m χ 0 1 , and cτ 0 (if applicable). The approach is based on the profile likelihood ratio as the test statistic [93], the CL s criterion [94,95], and the asymptotic formulae [96] to approximate the distributions of the test statistic under the SM-background-only and signal+background hypotheses. An Asimov data set [96] is used to determine the expected σ UL on the allowed cross section for a given model. Potential signal contributions to event counts in all bins of the SR and CRs are considered.
The experimental acceptance times efficiency (Aε) is evaluated independently for each model, defined in terms of m SUSY , m χ 0 1 , and cτ 0 (if applicable). The effects of several sources of uncertainty in Aε, as well as the potential for migration of events between bins of the SR, are considered. Correlations are taken into account where appropriate, including those relevant to signal contamination that may contribute to counts in the CRs.
The statistical uncertainty arising from the finite size of simulated samples can be as large as ≈30%. The Aε for models with a compressed mass spectrum relies on jets arising from ISR, the modelling of which is evaluated using the technique described in Ref. [67]. The associated uncertainty can be as large as ≈30%. The corrections to the jet energy scale (JES) evaluated with simulated events can lead to variations in event counts as large as ≈25% for models yielding high jet multiplicities. The uncertainties in the modelling of scale factors applied to simulated event samples that correct for differences in the b tagging efficiency and mistag probabilities can be as large as ≈20%. Table 4 defines a number of benchmark models that are close to the limit of the search sensitivity. All model families are represented, and the model parameters (m SUSY , m χ 0 1 , and cτ 0 if applicable) are chosen to select models with large and small differences in m SUSY and m χ 0 1 , as well as a range of cτ 0 values. Table 4 summarizes the aforementioned uncertainties for each benchmark model, presented in terms of a characteristic range that is representative of the variations observed across the bins of the SR. The upper bound for each range may be subject to moderate statistical fluctuations.
Additional subdominant contributions to the total uncertainty are also considered. The uncertainty in the integrated luminosity is determined to be 2.5% [23]. Uncertainties in the production cross section arising from the choice of the PDF set, and variations therein, as well as vari- Table 4: A list of benchmark simplified models organized according to production and decay modes (family), the Aε, representative values for some of the dominant sources of systematic uncertainty, and the expected and observed upper limits on the production cross section σ UL relative to the theoretical value σ th calculated at NLO+NLL accuracy. Additional uncertainties concerning the T1qqqqLL models are not listed here and are discussed in the text. ations in µ F and µ R at LO are considered. Uncertainties in event migration between bins from variations in the PDF sets are assumed to be correlated with, and adequately covered by, the uncertainties in the modelling of ISR. Uncertainties from µ F and µ R variations are determined to be ≈5%. The effect of a 5% uncertainty in the total inelastic cross section [91] is propagated through the weighting procedure that corrects for differences between the simulated and measured pileup, resulting in event count variations of ≈10%. Uncertainty in the modelling of the efficiency to identify high-quality, isolated leptons is ≈5% and is treated as anticorrelated between the SR and µ+jets and µµ+jets CRs. The uncertainty in the trigger efficiency to record signal events is <10%.
The Aε for the T1qqqqLL family of models depends on cτ 0 in addition to m g and m χ 0 1 . Scenarios involving a compressed mass spectrum or gluinos with cτ 0 10 m increase the probability that the decay of the gluino-pair system escapes detection, and the Aε is reduced for such models, as indicated in Table 4, because of an increased reliance on jets from ISR. Scenarios with m g − m χ 0 1 100 GeV and 1 cτ 0 10 m often lead to one or both gluinos decaying within the calorimeter systems to yield energetic jets comprising particle candidates that have no associated charged particle track. Hence, the efficiencies for the event vetoes related to the jet identification and f j 1 h ± requirements, described in Sections 2 and 3.1, can be as low as ≈90% and ≈30%, respectively, for this region of the model parameter space. Uncertainties as large as ≈10% are assumed. The efficiencies for all other scenarios are typically ≈100%. Jet identification requirements in the trigger logic lead to inefficiencies and uncertainties not larger than 2%. Finally, models with 1 cτ 0 10 mm often lead to jets that are tagged by the CSV algorithm with efficiencies as high as ≈60%, which are comparable to the values obtained for jets originating from b quarks. Uncertainties of 20-50% in the tagging efficiency are assumed to cover differences with respect to jets originating from b quarks, as indicated in Table 4. Figure 3 summarizes the excluded regions of the mass parameter space for the nine families of simplified models. The regions are determined by comparing σ UL with the theoretical cross sections σ th calculated at NLO+NLL accuracy. The former value is determined as a function of m SUSY and m χ 0 1 , while the latter has a dependence solely on m SUSY . The exclusion of models is evaluated using observed data counts in the signal region (solid contours) and also expected counts based on an Asimov data set (dashed contours). The observed excluded regions for the T1bbbb and T1tttt families, as shown in Fig. 3 (lower), can be up to 2-3 standard deviations weaker than the expected excluded regions when m g − m χ 0 1 ≈ 350 GeV. These differences are typically due to fluctuations in data for events that satisfy n jet ≥ 5 and n b ≥ 3. Figure 3 (lower) also allows a comparison of the sensitivity to T1qqqq and T1qqqqLL models, which assume the prompt-decay and metastable gluino scenarios, respectively. The latter scenario leads to a monojet-like final state as the gluino escapes detection, resulting in a reach in m g that is independent of m χ 0 1 . Figure 4 summarizes the evolution of the search sensitivity to the T1qqqqLL family of models as a function of cτ 0 . Each subfigure presents the observed σ UL as a function of m g and m χ 0 1 for simplified models that involve the production of gluino pairs. The excluded mass regions based on the observed and expected values of σ UL are also shown, along with contours determined under variations in theoretical and experimental uncertainties. The top row of subfigures cover the range 1 < cτ 0 < 100 µm and demonstrate coverage comparable to the T1qqqq promptdecay scenario. A moderate improvement in sensitivity for models with 1 cτ 0 10 mm is observed because of the additional signal-to-background discrimination provided by the n b variable. The sensitivity is reduced for models with lifetimes in the region cτ 0 > 100 mm because of a lower acceptance for the jets from the gluino decay and an increased reliance on jets from ISR. The coverage is independent of cτ 0 beyond values of 10 m and comparable to the limiting case of a metastable gluino.
A nonnegligible fraction of R-hadrons that traverse the muon chambers before decaying are identified as muons by the PF algorithm. The fraction is dependent on the R-hadron model and the choice of parameters that affect the hadronization model and matter interactions. The signal Aε is strongly dependent on cτ 0 due to the muon veto employed by this search. Under these assumptions, the excluded mass regions shown in Fig. 4 weaken by 50-200 GeV for models with cτ 0 1 m, with the largest change occurring at cτ 0 ≈ 10 m. The change is negligible for models with cτ 0 below 1 m. Table 5 summarizes the strongest expected and observed mass limits for each family of models. The simplified result based on the 32-bin schema, summarized in Appendix A, yields limits on σ UL that are typically a factor ≈2 weaker than those obtained with the nominal result.   Figure 3: Observed and expected mass exclusions at 95% CL (indicated, respectively, by solid and dashed contours) for various families of simplified models. The upper subfigure summarises the mass exclusions for five model families that involve the direct pair production of squarks. The first scenario considers the pair production and decay of bottom squarks (T2bb). Two scenarios involve the production and decay of top squark pairs (T2tt and T2cc). The grey shaded region denotes T2tt models that are not considered for interpretation. Two further scenarios involve, respectively, the production and decay of light-flavour squarks, with different assumptions on the mass degeneracy of the squarks as described in the text (T2qq 8fold and T2qq 1fold). The lower subfigure summarises three scenarios that involve the production and prompt decay of gluino pairs via virtual squarks (T1bbbb, T1tttt, and T1qqqq). A final scenario involves the production of gluinos that are assumed to be metastable on the detector scale (T1qqqqLL).

Summary
A search for supersymmetry with the CMS experiment is reported, based on a data sample of pp collisions collected in 2016 at √ s = 13 TeV that corresponds to an integrated luminosity of 35.9 ± 0.9 fb −1 . Final states with jets and significant missing transverse momentum p miss T , as expected from the production and decay of massive gluinos and squarks, are considered. Signal events are categorized according to the number of reconstructed jets, the number of jets identified as originating from bottom quarks, and the scalar and vector sums of the transverse momenta of jets. The standard model background is estimated from a binned likelihood fit to event yields in the signal region and data control samples. The observed yields in the signal region are found to be in agreement with the expected contributions from standard model processes. Supplemental material is provided to aid further interpretation of the result in Appendix A.
Limits are determined in the parameter spaces of simplified models that assume the production and prompt decay of gluino or squark pairs. The strongest exclusion bounds (95% confidence level) for squark masses are 1050, 1000, and 1325 GeV for bottom, top, and massdegenerate light-flavour squarks, respectively. The corresponding mass bounds on the neutralino χ 0 1 from squark decays are 500, 400, and 575 GeV. The gluino mass is probed up to 1900, 1650, and 1650 GeV when the gluino decays via virtual states of the aforementioned squarks. The strongest mass bound on the χ 0 1 from the gluino decay is 1150 GeV. Sensitivity to simplified models inspired by split supersymmetry is also demonstrated. These models assume the production of long-lived gluino pairs that decay to final states containing displaced jets and p miss T from the undetected χ 0 1 particles. The long-lived gluino, with an assumed proper decay length cτ 0 , is expected to hadronize with SM particles and form a bound state known as an R-hadron. The model-dependent matter interactions of R-hadrons are not considered by default. The sensitivity of this search is only moderately dependent on these matter interactions for models with cτ 0 1 m, while no dependence is found for models with cτ 0 below 1 m. Models that assume a χ 0 1 mass of 100 GeV and gluino masses up to 1600 GeV are excluded for a proper decay length cτ 0 below 0.1 mm. The bound on the gluino mass strengthens to 1750 GeV at cτ 0 = 1 mm, before weakening to 900-1000 GeV for models with cτ 0 > 10 m. For all values of cτ 0 considered, the exclusion bounds on the gluino mass weaken to about 1 TeV when the difference between the gluino and χ 0 1 mass is small. The search provides cov-erage of the cτ 0 parameter space for models involving long-lived gluinos, such as the region cτ 0 1 mm, that is complementary to the coverage provided by dedicated techniques at the LHC.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centres and personnel of the Worldwide LHC Computing Grid for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC and the CMS detector provided by the following funding agencies: BMWFW and FWF (Aus     [46] ATLAS Collaboration, "Search for long-lived, massive particles in events with displaced vertices and missing transverse momentum in √ s = 13 TeV pp collisions with the ATLAS detector", (2017). arXiv:1710.04901. Submitted to Phys. Rev. D.
[52] CMS Collaboration, "Performance of CMS muon reconstruction in pp collision events at √ s = 7 TeV", JINST 7 (2012) P10002, doi:10.1088/1748-0221/7/10/P10002, arXiv:1206.4071. Table 6: Summary of the nominal (n jet , n b , H T , H miss T ) binning schema. Each entry (and the following entry, if present) signifies the lower (upper) bound of an H miss T bin within a given (n jet , n b , H T ) bin. Unique or final entries represent H miss T bins unbounded from above. A dash (-) signifies that the H T bin in a given (n jet , n b ) category is not used in the analysis, in which case counts in high-H T bins are integrated into the adjacent lower-H T bin. For monojet events, H T ≡ H miss T . The a denotes asymmetric p T thresholds for the two highest p T jets.