Search for long-lived particles using displaced jets in proton-proton collisions at $\sqrt{s} = $ 13 TeV

An inclusive search is presented for long-lived particles using displaced jets. The search uses a data sample collected with the CMS detector at the CERN LHC in 2017 and 2018, from proton-proton collisions at a center-of-mass energy of 13 TeV. The results of this search are combined with those of a previous search using a data sample collected with the CMS detector in 2016, yielding a total integrated luminosity of 132 fb$^{-1}$. The analysis searches for the distinctive topology of displaced tracks and displaced vertices associated with a dijet system. For a simplified model, where pair-produced long-lived neutral particles decay into quark-antiquark pairs, pair production cross sections larger than 0.07 fb are excluded at 95% confidence level (CL) for long-lived particle masses larger than 500 GeV and mean proper decay lengths between 2 and 250 mm. For a model where the standard model-like Higgs boson decays to two long-lived scalar particles that each decays to a quark-antiquark pair, branching fractions larger than 1% are excluded at 95% CL for mean proper decay lengths between 1 mm and 340 mm. A group of supersymmetric models with pair-produced long-lived gluinos or top squarks decaying into various final-state topologies containing displaced jets is also tested. Gluino masses up to 2500 GeV and top squark masses up to 1600 GeV are excluded at 95% CL for mean proper decay lengths between 3 and 300 mm. The highest lower bounds on mass reach 2600 GeV for long-lived gluinos and 1800 GeV for long-lived top squarks. These are the most stringent limits to date on these models.

Given the large variety of the BSM scenarios that lead to displaced-jets signatures, it is important to make the displaced-jets search as model independent as possible. In this paper, we present an inclusive search for LLPs decaying into jets, with at least one LLP having a decay vertex within the tracker acceptance, but which is displaced from the production vertex by up to 550 mm in the plane transverse to the beam direction. The search looks for a pair of jets known as dijets, where the jets are clustered from energy deposits in the calorimeters. For jets arising from the decay of an LLP, the associated tracks are usually displaced from the primary vertices (PVs), and the decay vertex can be reconstructed from the displaced tracks. The properties of the tracks and the decay vertex can provide discrimination power to distinguish long-lived signals from SM backgrounds. As mentioned above, a large number of models predict LLPs decaying into displaced jets. Our tests for some of these will be discussed in detail in Section 3.
Events used in this analysis were collected with the CMS detector [39] at the LHC from protonproton (pp) collisions at a center-of-mass energy of 13 TeV in 2017 and 2018, corresponding to an integrated luminosity of 95.9 fb −1 . The results are combined with those of a previous displaced-jets search using the events collected in 2016 [40], yielding a total integrated luminosity of 132 fb −1 . For the models that were not studied in the 2016 displaced-jets search, additional simulated signal samples have been produced following the 2016 run condition of the CMS detector. These additional samples are then processed with the reconstruction and selection procedures described in Ref.
[40] to compute the additional signal yields and systematic uncertainties for the 2016 data that are used in the combination.
Compared to the 2016 displaced-jets search, a set of new techniques that significantly improves the sensitivity to long-lived signatures is implemented in this analysis. The new techniques include one additional dedicated trigger aimed at selecting events containing displaced jets to recover efficiencies for high-mass LLPs, an auxiliary nuclear interactions (NIs) veto map to improve background rejection, a dedicated variable based on the sum of signed impact parameters of the tracks assigned to the displaced vertex, and the use of machine learning techniques to improve signal-to-background discrimination. With these new techniques, compared to the 2016 search, we have reduced the background rate by approximately a factor of three, while significantly increasing the signal efficiencies for almost all signal points in different LLP models. Results of searches for similar LLP signatures with hadronic decays at √ s = 13 TeV have also been reported by ATLAS [41][42][43][44][45] and CMS [46][47][48].
The paper is organized as follows. A brief description of the CMS detector is introduced in Section 2. The data and the simulated samples are described in Section 3. Section 4 details the event reconstruction and the preselection criteria. Section 5 describes the event selections and the background estimation methods. The systematic uncertainties are summarized in Section 6. The observation and the interpretation of the results are described in Section 7. The paper is summarized in Section 8.

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid.
The silicon tracker measures charged particles within the pseudorapidity range |η| < 2.5. During the LHC run in 2017 and 2018, the silicon tracker consisted of 1856 silicon pixel and 15 148 silicon strip detector modules, and it occupies a cylindrical volume around the interaction point with a length of 5.8 m and a diameter of 2.6 m. For nonisolated particles with 1 < p T < 10 GeV and |η| < 1.4, the track resolutions are typically 1.5% in p T , and 25-75 µm in the transverse impact parameter [49].
In the region |η| < 1.74, the HCAL cells have widths of ∆η = 0.087 in pseudorapidity and ∆φ = 0.087 in azimuth. In the η-φ plane, and for |η| < 1.48, the HCAL cells map on to 5×5 arrays of ECAL crystals to form calorimeter towers projecting radially outward from the nominal interaction point. For |η| > 1. 74, the coverage of the towers increases progressively to a maximum of 0.174 in ∆η and ∆φ. Within each tower, the energy deposits in ECAL and HCAL cells are summed to define the calorimeter tower energies, and are subsequently used to provide the energies and directions of hadronic jets.
Events of interest are selected using a two-tiered trigger system [50]. 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 (HLT), 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 around 1 kHz before data storage.
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [39].

Datasets and simulated samples
Data were collected with two dedicated triggers aimed at selecting events containing displaced jets. At the HLT, jets are reconstructed from the energy deposits in the calorimeter towers, clustered using the anti-k T algorithm [51,52] with a distance parameter of 0.4. In this process, the contribution from each calorimeter tower is assigned a momentum, the magnitude and the direction of which are given by the energy measured in the tower and the coordinates of the tower. The raw jet energy is obtained from the sum of the tower energies, and the raw jet momentum from the vector sum of the tower momenta, which results in a nonzero jet mass. The raw jet energies are then corrected [53] to establish a uniform relative response of the calorimeter in η and a calibrated absolute response in transverse momentum p T .
Identification of the PV is a prerequisite for the selection of displaced jets at the HLT. Events may contain multiple PVs, corresponding to multiple pp collisions occurring in the same bunch crossing. The candidate vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm [51,52] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets. More details are given in Section 9.4.1 of Ref. [54].
The first trigger, referred to as the "displaced" trigger, requires H T > 430 GeV, where H T is the scalar sum of the jet p T for all jets that have p T > 40 GeV and |η| < 2.5 in the event. The trigger also requires the presence of at least two jets, with the following requirements satisfied for each jet: • p T > 40 GeV and |η| < 2.0; • at most two associated prompt tracks with p T > 1 GeV, where prompt tracks are those having a transverse impact parameter (IP 2D ) with respect to the leading PV smaller than 1.0 mm; and • at least one associated displaced track with p T > 1 GeV, where a displaced track is a track having an IP 2D larger than 0.5 mm and an impact parameter significance (Sig[IP 2D ]) larger than 5.0. The significance is defined as the ratio between the impact parameter and its uncertainty.
The second trigger, referred to as the "inclusive" trigger, requires H T > 650 GeV and the presence of at least two jets, each of them satisfying: • p T > 60 GeV and |η| < 2.0; and • at most two associated prompt tracks with p T > 1 GeV.
The displaced trigger is more efficient for low-mass LLPs, while the inclusive trigger is designed to recover the trigger efficiency for high-mass LLPs with small ( 3 mm) or large ( 300 mm) mean proper decay lengths (cτ 0 ).
The background sources in this search include NIs between outgoing particles and detector material, long-lived SM hadrons, and misreconstructed displaced vertices formed by accidentally crossing tracks. The background events mainly arise from SM events containing jets produced through the strong interaction, referred to as quantum chromodynamics (QCD) multijet events. The QCD multijet Monte Carlo (MC) sample is simulated at leading order with MADGRAPH5 aMC@NLO 2.4.2 [55]. Parton showering and hadronization are simulated with PYTHIA 8.226 [56]. The matching of jets from the matrix element calculations and parton shower jets is achieved using the MLM algorithm [57]. The PYTHIA parameters for the underlying event modeling are set to be the CP5 tune [58]. The set of parton distribution functions (PDFs) used for the production is the NNPDF3.1 NNLO PDF set [59]. The QCD multijet MC sample is mainly used to inspire the analysis strategy and to estimate systematic uncertainties, while the background estimation for this search is purely determined from data.
Feynman diagrams for the benchmark models studied in this paper are summarized in Fig. 1.
One of the benchmark signal models is a simplified model, where long-lived scalar neutral particles X are pair produced through a scattering process mediated by an off-shell Z boson. In this model, each X particle decays to a quark-antiquark pair, assuming equal branching fractions to u, d, s, c, and b quark pairs. The decays to top quark pairs are excluded to provide a simple final-state topology for this model, but it was important that the analysis strategy would still be sensitive to a variety of other models. We checked the impact on the signal efficiencies of excluding decays to top quark pairs, and found it to be small, where the relative changes of the signal efficiencies are generally at the order of a few percents. This model is referred to as the jet-jet model. The samples are produced with different resonance masses ranging from 50 to 1500 GeV, and with different proper decay lengths ranging from 1 to 10 4 mm.
Another signature we consider is the case where LLPs arise from exotic decays of an SM-like Higgs boson, which can happen in many BSM scenarios (a review can be found in the Section IV.6.6 of Ref. [60]), including "hidden Valley" models [18,19], twin Higgs models [36], and the folded SUSY model [61]. For the simulation, we use POWHEG 2.0 [62][63][64][65] to generate events containing a 125 GeV Higgs boson produced through gluon-gluon fusion. The 125 GeV Higgs boson then decays to two long-lived scalar particles S, and each scalar particle then decays to a quark-antiquark pair. Two scenarios are considered; in the first scenario the scalar particle has a branching fraction of 100% to decay to a down quark-antiquark pair, while in the second one the scalar particle has a branching fraction of 100% to decay to a bottom quark-antiquark pair.
The samples are produced with the scalar particle mass m S set to be 15, 40, or 55 GeV, while the cτ 0 of S varies from 1 to 3000 mm.
We also consider a group of SUSY models with different final-state topologies. The first one is a GMSB SUSY model [66] in the general gauge mediation scenario [14,15], where gluinos are pair produced and the gravitino is the lightest SUSY particle, while the gluino is the nextto-lightest supersymmetric particle. After the gluino is produced, it decays to a gluon and a gravitino, producing a single displaced jet and missing transverse momentum. This decay is suppressed by the SUSY-breaking scale, and therefore the gluino is long lived. The model in which this process occurs is referred to as the g → g G model. The samples are produced with gluino masses from 800 to 2500 GeV, and with the cτ 0 of the gluino varying from 1 to 10 4 mm.
The second SUSY model we consider is a mini-split SUSY model [6,7], referred to as the g → qq χ 0 1 model. In this model the gluino decays to a quark-antiquark pair and the lightest neutralino ( χ 0 1 ), with equal branching fractions to u, d, s, and c quark pairs. This decay is mediated by a squark, which is much heavier than the gluino. The squark's large mass suppresses the gluino decay, making it long lived. The mass of the neutralino is assumed to be 100 GeV, the samples are produced with gluino masses from 1400 to 3000 GeV, and the cτ 0 of the gluino varies from 1 to 10 4 mm. The third SUSY model is an RPV SUSY model [67] with minimal flavor violation, where gluinos are pair produced and long lived. Each long-lived gluino decays to top, bottom, and strange antiquarks through the RPV coupling λ 323 and the mediation of a virtual top squark [12], leading to a multijet final-state topology. This model is referred to as the g → tbs model. The samples are produced with gluino masses from 1200 to 3000 GeV, and a cτ 0 varying from 1 to 10 4 mm.
We also consider two other RPV SUSY models [68] with semileptonic decays, in which longlived top squarks are pair produced, and each top squark decays to a bottom quark (down quark) and a charged lepton via RPV couplings λ 133 , λ 233 , and λ 333 (λ 131 , λ 231 , and λ 331 ) [12]. The decay rate to each of the three lepton flavors is assumed to be equal. The two models are referred to as the t → b ( t → d ) models. The samples are produced with different top squark masses from 600 to 2000 GeV, and a cτ 0 varying from 1 to 10 4 mm.
Finally, we consider another SUSY model, referred to as the t → dd model, motivated by dynamical RPV (dRPV) [69,70], where long-lived top squarks are pair produced, and each top squark decays to two down antiquarks via a nonholomorphic RPV coupling η 311 [71]. The nonholomorphic RPV coupling is suppressed by some large scale M, thus giving rise to long  Figure 1: The Feynman diagrams for the different long-lived models considered, including the jet-jet model (upper left), models with an exotic decay of the SM-like Higgs boson (upper right), general gauge mediation models with g → g G decay (second row, left), mini-split SUSY with g → qq χ 0 1 decay (second row, right), RPV SUSY with g → tbs decay (third row, left), RPV SUSY with t → b decay (third row, right), RPV SUSY with t → d decay (lower left), and dRPV SUSY with t → dd decay (lower right).
lifetimes. The samples are produced with different top squark masses from 800 to 1800 GeV, and a cτ 0 varying from 1 to 10 4 mm. PYTHIA 8.226 is used for the production of the signal samples, and the PDF set used for the production is NNPDF3.1 LO. For SUSY-particle production, the PYTHIA 8.226 simulation is cross checked with MADGRAPH5 aMC@NLO 2.4.2 for representative signal points, where the MAD-GRAPH5 aMC@NLO simulation is performed at LO with up to two additional outgoing partons. The resulting signal efficiencies are found to be consistent within the statistical uncertainties. The PYTHIA parameters for the underlying event modeling are set to be the CP2 tune [58]. In the SUSY models, a long-lived gluino or top squark can form a hadronic state through strong interactions, an R-hadron [9,72,73], which is simulated with PYTHIA. The interactions of the R-hadron with matter were studied following the simulation described in Ref. [74,75], and were found to have negligible impact on this analysis, since they have very little influence on the vertex reconstruction.
The simulated background and signal events are processed with a GEANT4-based [76] simulation for the detailed CMS detector response. To take into account the effects of additional pp interactions within the same or nearby bunch crossings ("pileup"), additional minimum-bias events are overlaid on the simulated events to match the pileup distribution observed in the data.

Event reconstruction and preselection
This search examines dijet candidates in a given event. The algorithms for the offline jet reconstruction and PV selection are the same as those applied at the HLT (as described in Section 3), except that the full offline information is used. To make sure that the online H T and jet p T requirements in the displaced-jet triggers reach full efficiency, we apply selections on the offline H T of the event as well as on the p T and η of each jet. After the trigger selection, if an event passes the displaced trigger, we require the event to have offline H T > 500 GeV, and dijet candidates are formed from all possible pairs of jets in the event, with the jets satisfying p T > 50 GeV and |η| < 2.0. On the other hand, if an event only passes the inclusive trigger, it is required to have offline H T > 700 GeV, and the dijet candidates are formed from all possible pairs of jets in the event, with the jets satisfying p T > 80 GeV and |η| < 2.0.
In this search, the track candidates are required to have p T > 1 GeV and to be high-purity tracks. The high-purity selection is based on track information (such as the normalized χ 2 of the track fit, the impact parameters, and the number of hits in different tracker layers) to reduce the fraction of misreconstructed tracks, and the selection is optimized separately for each iteration of the tracking [49], so that it is efficient for selecting tracks with different displacements. More details of the high-purity selection can be found in Ref. [49]. The η and φ of a given track are determined by the direction of its momentum vector at the closest approach point to the leading PV. For a given dijet candidate, we associate track candidates with each jet by requiring that ∆R < 0.5, where ∆R = √ (∆η) 2 + (∆φ) 2 is the angular distance between the jet axis and the track direction. When a track satisfies ∆R < 0.5 for both jets, it is associated with the jet giving the smaller ∆R.
After associating track candidates with each jet, the next step is to reconstruct a secondary vertex (SV) for each dijet candidate. From all the tracks associated with a dijet candidate, we select displaced tracks that satisfy IP 2D > 0.5 mm and Sig[IP 2D ] > 5.0. We then attempt to reconstruct an SV from these displaced tracks using an adaptive vertex fitter algorithm [77]. The reconstructed SV is not required to have associated tracks from both jets, so that the search can be sensitive to the models where the LLP decays to a single displaced jet. To improve the signal-to-background discrimination, we implement a set of preselection criteria on the dijet/SV candidates, which are described in the rest of this section.
To ensure the quality of the vertex reconstruction, the SV is selected only if it is reconstructed with a χ 2 per degree of freedom (χ 2 /n dof ) of less than 5.0. In order to suppress long-lived SM mesons and baryons, the invariant mass of the vertex is required to be larger than 4 GeV, and the transverse momentum of the vertex is required to be larger than 8 GeV, where the fourmomentum of the vertex is calculated assuming the charged pion mass for all assigned tracks.
We only consider dijet candidates that have a reconstructed SV satisfying the above requirements. Furthermore, SVs in background events tend to have only one track with a high value for IP 2D , corresponding to the tail of the impact parameter distribution. We therefore consider the track with the second-highest Sig[IP 2D ] among the tracks that are assigned to the SV, since this provides a more sensitive discriminant for identifying displaced jets. We require the second-highest Sig[IP 2D ] to be larger than 15.
We also compute another quantity , which is the ratio between the sum of energy for all the tracks assigned to the SV and the sum of the energy for all the tracks associated with the two jets: Since is expected to be large for displaced-jet signatures, dijet candidates with smaller than 0.15 are rejected.
An additional variable, ζ, is defined to characterize the contribution of prompt activity to the jets. For each track associated with a jet, we identify the PV (including the leading PV and the pileup vertices) with the minimum three-dimensional (3D) impact parameter significance to the track. If this minimum value is smaller than 5, we assign the track to this PV. Then for each jet, we compute the track energy contribution from each PV, and the PV with the largest track energy contribution to the jet is chosen. Finally, we define ζ as where ∑ track∈PV i E Jet i track is the sum of the track energy coming from the most compatible PV for a given jet, while E Jet i is the energy of a given jet, thus ζ is the charged energy fraction of the dijet associated with the most compatible PVs. For displaced-jet signatures, ζ tends to be small since the jets are not compatible with PVs. Dijet candidates with ζ larger than 0.2 are rejected.
To suppress the background events arising from NIs in the tracker material, we compare the positions of the SVs with a map of the distribution of material in the inner tracker. The map was obtained from the distribution of NI candidate vertices, which are reconstructed using the adaptive vertex fitter on a sample of events collected with isolated single-muon triggers. The NI candidates are required to satisfy the following criteria: • The tracks are required to be associated with dijet candidates, which are formed from the jets having p T > 10 GeV and |η| < 2.0; • The associated tracks must have p T > 0.2 GeV, high purity, IP 2D > 0.5 mm, and Sig[IP 2D ] > 5.0; • vertex track multiplicity is larger than 3; • the ratio of the energy sum for SV tracks to that for all tracks is less than 0.15 for the dijet candidate; • vertex L xy significance is larger than 200, where the transverse decay length L xy is the distance between the SV and the leading PV; and • vertex χ 2 /n dof is smaller than 3.0.
After these selections, the distribution of the NI vertex candidates is transferred to an NI-veto map in the transverse plane, with |x| and |y| < 25 cm, as shown in Fig. 2. To suppress misreconstructed vertices and the displaced vertices produced by decaying long-lived SM mesons and baryons, we only select the region where the NI vertex density is above a threshold that varies for different layers of the pixel detector. In the NI-veto map, we can clearly see the structures of the beam pipe (at r = √ x 2 + y 2 ≈ 23 mm), the four pixel layers (at r ≈ 29, ≈ 68, ≈ 109, and ≈ 160 mm), and the support rails (at r ≈ 200 mm). In our search, any SV candidate that overlaps with the NI-veto map is rejected. The loss of the fiducial volume within r < 300 mm due to the veto is around 4%, and the efficiencies for signal events to pass this selection are generally well above 90%. In the veto no requirement is placed on the z coordinates of the SVs, but the impact of restricting the veto to the barrel region of the pixel detector (|z| < 27 cm) is negligible on the signal efficiencies. A similar study on the structure of the CMS inner tracking system using a more sophisticated NI reconstruction technique with 2016 data has been reported in Ref. [78].  [79]. The structures of the different pixel layers can be clearly seen. Right: the efficiency for a given vertex candidate to pass the NI-veto as a function of radius r.
The preselection criteria for this search, summarized in Table 1, are efficient for a wide range of long-lived models with different final-state topologies.

Event selection and background prediction
After reconstructing the SV using the adaptive vertex fitter, we employ an auxiliary algorithm to check the consistency between the SV system and the dijet system. For each displaced track   [80]. During the clustering, two clusters are merged when the smallest L exp xy difference between the two clusters is smaller than 15% of the L xy of the SV. After the clustering procedure is finished, if more than one cluster is formed, the one closest to the SV is selected. The cluster root-mean-square (RMS), taken to be the relative RMS of individual tracks L exp xy with respect to the SV L xy , is computed to provide signal-to-background discrimination: where N tracks is the number of tracks in the selected cluster.
For each track assigned to the SV, a sign is given to the IP 2D and Sig[IP 2D ] based on the angle between the dijet direction and the impact parameter vector that points from the leading PV to the closest approach point (with respect to the leading PV) of the track in the transverse plane. The sign is positive if this angle is smaller than π/2; otherwise the sign is negative. A new variable, κ, is then introduced as the signed Sig[IP 2D ] sum of the six leading tracks from the SV (where the tracks are ordered by the absolute values of their Sig[IP 2D ]): If the track multiplicity is smaller than six, the sum is taken over all the tracks from the SV. For background processes, since the tracks assigned to the SV are uncorrelated with the dijet direction, the signed Sig[IP 2D ] of different tracks tend to cancel each other, therefore κ peaks sharply around zero. On the other hand, for displaced jets that originate from the SV, the directions of the tracks will be highly correlated with the dijet direction, therefore κ is significantly different from zero and |κ| tends to be large.
To improve signal-to-background discrimination and to define a region with signal events enriched, we proceed to construct a multivariate discriminant based on the following variables for the vertex/dijet candidates: • Vertex track multiplicity; • Vertex L xy significance; • Cluster RMS; • The magnitude of the signed Sig[IP 2D ] sum of the six leading tracks |κ|.
The distributions of the four variables are shown in Fig. 3, with displaced-jet triggers, offline H T , and offline jet kinematic variables selections (described in Section 4) applied. For the multivariate discriminant we utilize the gradient boosted decision tree (GBDT) algorithm [81][82][83], with cross entropy as the loss function. The GBDT algorithm is implemented using the TMVA (toolkit for multivariate data analysis) package [84] interfaced with SCIKIT-LEARN [85]. Given the large cross section of the QCD multijet process and the relatively low H T threshold of our displaced-jet triggers, the event count of the simulated QCD multijet sample (after preselections) is insufficient for the GBDT training, since it is much smaller than the number of expected QCD multijet events in the analyzed data sample. Therefore, for the background sample in the GBDT training, we use the data in the following region: • events are selected by the displaced-jet triggers, and pass the offline H T and jet kinematic variables selections; • < 0.12 for the dijet candidate, making this region orthogonal to the signal region; • the veto using the NI-veto map is not applied; • all the other preselection criteria are satisfied.
For the signal sample in the GBDT training, simulated jet-jet model events that pass the preselection criteria are used, with m X = 100, 300, and 1000 GeV, and with cτ 0 = 1, 10, 100, 1000 mm. If there is more than one dijet/SV candidate passing the selection criteria in a given event, the one with the largest track multiplicity is chosen for the training. If the track multiplicities are the same, the one with the smallest SV χ 2 /n dof is chosen. An event weight is assigned separately to each signal point with a given m X and cτ 0 such that the sum of weights is identical for each point, thus each signal point has the same priority in the training. Twenty percent of the events in the signal and background samples are randomly selected to validate the performance of the GBDT and to make sure it is not overtrained.
The GBDT output values or scores for data, simulated QCD multijet events, and simulated signal events are shown in Fig. 4. The signal efficiencies for this search are measured with simulated signal events produced separately. The background prediction is purely based on some other control samples in data, which are different from the one used for the GBDT training.
Although only the jet-jet model is used as the signal sample in the GBDT training, the GBDT is highly efficient in selecting signatures of other LLP models with different final-state topologies, since we have explicitly chosen the input variables to make the GBDT as model-independent as possible.
In addition to the GBDT score g, we use another variable N 3D tracks in the final event selection, which is the number of 3D prompt tracks in a single jet, where the 3D prompt tracks are the tracks that have 3D impact parameters with respect to the leading PV smaller than 0.3 mm.
If more than one dijet candidate passes the preselection criteria described in Section 4, the one with the largest GBDT score is selected. If the GBDT scores are the same, the one with the smallest SV χ 2 /n dof is selected. In the final signal region, the candidate is further required to pass three final selection criteria, which are: • Selection 1: for the leading jet in p T , N 3D tracks is smaller than 3; • Selection 2: for the subleading jet, N 3D tracks is smaller than 3; and ] sum |κ| (lower right), for data, simulated QCD multijet events, and simulated signal events. Data and simulated events are selected with the displaced-jet triggers and with the offline H T , jets p T , and η selections applied. For a given event, if there is more than one SV candidate being reconstructed, the one with the largest vertex track multiplicity is chosen. If the track multiplicities are the same, the one with the smallest χ 2 /n dof is chosen. The lower panels show the ratios between the data and the simulated QCD multijet events. The blue shaded error bands and vertical bars represent the statistical uncertainties. Three benchmark signal distributions are shown (dashed lines) for the jet-jet model with m X = 300 GeV and varying cτ 0 . For visualization purposes, each signal process is given a cross section that yields 10 6 events produced in the analyzed data sample. The distributions of the GBDT output score for data, simulated QCD multijet events, and simulated signal events. Data and simulated events are selected with the displaced-jet triggers and with the offline H T , jets p T , and η selections applied. For a given event, if there is more than one SV candidate being reconstructed, the one with the largest vertex track multiplicity is chosen. If the track multiplicities are the same, the one with the smallest χ 2 /n dof is chosen. The lower panel shows the ratio between the data and the simulated QCD multijet events. The blue shaded error bands and vertical bars represent the statistical uncertainties. Three benchmark signal distributions are shown (dashed lines) for the jet-jet model with m X = 300 GeV and varying cτ 0 . For visualization purposes, each signal process is given a cross section corresponding to 10 6 events produced in the analyzed data sample. The signal events shown in this plot are not used in the GBDT training.
• Selection 3: the GBDT score g is larger than 0.988.
The numerical values for the selection criteria are chosen by optimizing the discovery potential of 5 standard deviations based on the Punzi significance [86] for the jet-jet model and the g → g G model across different LLP masses and lifetimes. The chosen models encompass the displaced-dijet and displaced-single-jet signatures, with m X = 100, 300, and 1000 GeV for the jet-jet model, and with m g = 600, 1000, and 1600 GeV for the g → g G model, while cτ 0 is taken to be 1, 10, 100, and 1000 mm. Thus there are 24 signal points considered in total for the selection optimization.
Based on the three selection criteria, we can define eight nonoverlapping regions A-H, which include the final signal region. The event counts in different regions are N f f f , N p f f , · · · , and N ppp , for regions A, B, · · · , and H, respectively, as shown in Table 2. The region H is the region where the events pass all the three selection criteria, and thus is the final signal region. Events in the remaining regions (A-G) fail one or more of the three selection criteria. Since the three selection criteria have little correlation between them for the background events, a property that has been verified with simulated QCD multijet events, the background yield in the signal region H can be estimated by different ratios of event counts in regions A-G, where the ratio uses the fraction of events passing to those failing the GBDT selection (Selection 3) and is taken as the central value of the predicted background yields. Three additional ratios are computed using the events failing either one of the first two selections: Fail These cross-checks provide an important test of the robustness of the background prediction and the assumption that the three selection criteria are minimally correlated. Differences between the predictions obtained with the nominal ratio b nominal and the cross-checks are used to estimate the systematic uncertainties in the background prediction.
The method described above can be generalized to predict the background yield in arbitrary intervals of the GBDT score g, where the first two selections are also satisfied. In this way we can verify our background prediction method in the different bins of GBDT score g. The background prediction method is first tested with simulated QCD multijet events and simulated signal events, and is found to be robust both with and without signal contaminations. We then test the background prediction in data, for which we define a control region by inverting the selection on , requiring to be smaller than 0.15. In order to improve the statistical precision, we remove the selection requirement that uses the NI-veto map, although the background pre-diction method is still robust with the presence of the NI-veto map. The predicted background yields and observed events in the control region are shown in Fig. 5, which shows a good agreement between the predicted and observed yields. The predicted background yields in the different bins of GBDT score are correlated, since the events that are used for background predictions in lower bins are also used in the background predictions in higher bins.

Systematic uncertainties
The systematic uncertainty in the background prediction is taken to be the largest deviation of the three cross checks from the nominal prediction b nominal , and is found to be 52% in the final signal region where the GBDT score g is larger than 0.988.
The systematic uncertainties in the integrated luminosity for 13 TeV pp collision data are 2.3% and 2.5% for the 2017 and 2018 [87,88] data taking periods, respectively, and are modeled as uncorrelated nuisance parameters between the years.
The systematic uncertainty in the signal yields arising from the pileup modeling is estimated by varying the inelastic pp cross section [89] used in the pileup distribution determination by 4.6%. The resulting variations of the signal yields are generally smaller than 1% and are well within the statistical fluctuations, therefore the impact of the pileup modeling is negligible and no corresponding systematic uncertainty is assigned.
The systematic uncertainty in the signal yields due to the online H T requirements of the displacedjet triggers is estimated by comparing the efficiency of the online H T requirements measured in data with the one measured in the QCD multijet MC sample. The efficiencies are measured using the events collected with an isolated single-muon trigger. The discrepancies between the measurements in data and MC simulation are parameterized as functions of offline H T , and corresponding corrections are applied to the simulated signal events. The signal yields are then recalculated for different masses and mean proper decay lengths. The largest correction in the signal yield for a given model is taken as the corresponding systematic uncertainty, and is found to be smaller than 2%.
The systematic uncertainty in the signal yields due to the online jet p T requirements of the displaced-jet triggers is estimated by measuring the per-jet efficiencies of the online jet p T requirements in data and in the QCD multijet MC sample, using events collected with a prescaled H T trigger that requires H T > 425 GeV. Corrections for the discrepancies between the measurements in data and MC simulation are applied to the simulated signal events, and the signal yields are recalculated. The correction is more significant for low-mass LLPs (with m ∼ 50 GeV), and is smaller than 8%, which is taken as the corresponding systematic uncertainty.
To estimate the systematic uncertainty due to the online tracking requirements of the displacedjet triggers, we measure the per-jet efficiencies of the online tracking requirements as functions of the number of prompt tracks and displaced tracks, using events collected with the prescaled H T trigger. The efficiencies obtained in data are found to be consistent with the efficiencies obtained in MC simulations. Therefore no corresponding systematic uncertainty is assigned.
To estimate the impact of the possible mismodeling of the GBDT score in MC simulation on the signal yields, we compare the distribution of the GBDT score in simulated QCD multijet events with the distribution measured in data, using events collected with the prescaled H T trigger. The discrepancies between data and MC simulation are taken into account by varying the GBDT scores in the simulated samples by the same magnitude. The largest variation in the signal efficiency for a given model is taken as the corresponding systematic uncertainty in signal yields, and is found to be in the range of 4%-15%.
Similarly, the uncertainty in the signal yields due to impact parameter modeling is estimated by comparing the distribution of the impact parameters in simulated QCD multijet events with the distribution measured in data, also using events collected with the prescaled H T trigger. The impact parameters in simulated QCD multijet events are varied until the discrepancies between data and MC simulation are covered by the variation. The impact parameters in the simulated signal samples are then varied by the same magnitude. The largest variation in the signal efficiency for a given model is taken as the corresponding systematic uncertainty, and is found to be in the range of 8%-18%.
The impact of the jet energy scale uncertainty on the signal yields is estimated by varying the jet energy and p T by one standard deviation of the jet energy scale uncertainty [53]. The variations of the signal efficiencies are smaller than 3%, which is taken as the corresponding systematic uncertainty.
The impact of the PDF uncertainty on the signal acceptance is estimated by reweighting the simulated signal events with NNPDF, CT14 [90] and MMHT14 [91] PDF sets, and their associated uncertainty sets [92,93], following the PDF4LHC recommendation [92]. The uncertainty in signal efficiency for a given signal model is quantified by comparing the efficiencies calculated with alternative PDF sets and the ones with the nominal NNPDF set, and is found to be in the range of 4-6%.
The uncertainty in the signal yields due to the selection of the PV is estimated by replacing the leading PV with the subleading PV when calculating impact parameters and vertex displacement. The largest variation of the signal efficiency for a given signal model is taken as the corresponding systematic uncertainty, and is found to be in the range of 8-15%.
The various systematic uncertainties in the signal yields are summarized in Table 3.

Data in the signal region
The predicted background yields and the numbers of observed events in different GBDT ranges, after all the preselection criteria are applied, are shown in Fig. 6, with N 3D tracks smaller than 3 for both jets. The final signal region is defined by a GBDT score larger than 0.988, and the predicted background yield is 0.75 ± 0.44 (stat) ± 0.39 (syst) events in the data samples collected in 2017 and 2018. After all the preselection criteria are applied, the efficiencies for signal events to have at least one SV/dijet candidate satisfying GBDT > 0.988 are generally above 80-90% for the different LLP models considered with m LLP 300 GeV and 1 cτ 0 1000 mm, while the background rate has been reduced by a factor of around 3 × 10 3 by the GBDT selection. Event yields in data after different selection requirements have been applied are shown in Table 4. We observe one event in the final signal region, which is consistent with the predicted background yield. The observed event has an offline H T of 570 GeV, and an SV candidate with an L xy of ≈26 cm and 8 associated tracks. The position of the SV is close to one of the silicon strip layers, and was likely produced by NIs with the silicon strip detector.

Interpretation of the results
The signal yields in the final signal region are used to set limits on a variety of models. The signal efficiencies for representative signal points in different models can be found in Tables A.1 Figure 6: The predicted background yields and the number of observed events for the data in the signal region, with N 3D tracks smaller than 3 for both jets, shown for different bins of the GBDT scores. The background predictions in different bins are correlated, since the events that are used for background predictions in lower bins are also used in the background predictions in higher bins. For comparison, three benchmark signal points are also shown (dashed lines) for the jet-jet model with m X = 300 GeV and different lifetimes. For visualization purposes, each signal process is given a cross section that yields 100 events produced in the analyzed data sample. not studied in the 2016 displaced-jets search, we have produced additional signal MC samples following the 2016 run condition of the CMS detector. We then process these samples with the reconstruction and selection procedures implemented in the 2016 search to compute the additional signal yields and systematic uncertainties for the 2016 data that are used in the combination.
Upper limits on the cross section for a given model are presented with different masses and lifetimes by computing the 95% confidence level (CL) associated with each signal point using the CL s prescription [94][95][96][97], for which an LHC-style profile likelihood ratio [96,97] is taken as the test statistic. Systematic uncertainties are incorporated through the use of nuisance parameters treated according to the frequentist paradigm. The asymptotic approximation [96] is used for the calculation of the CL s values, which have been verified with full-frequentist results for representative signal points. Since the background yields of this search are small, the impact of the associated statistical or systematic uncertainties on the upper limits are also small.
The expected and observed upper limits on the pair production cross section for the jet-jet model at different values of cτ 0 and LLP mass m X are shown in Fig. 7, where a branching fraction of 100% for X to decay into a quark-antiquark pair is assumed. For a fixed LLP mass m X , the limits are most restrictive for cτ 0 between 3 and 300 mm. For smaller cτ 0 , the limits become less stringent, because we only select displaced tracks to reconstruct SVs and we veto dijet candidates with a large number of prompt tracks. The limits also become less restrictive for larger cτ 0 (cτ 0 > 300 mm), because the tracking efficiency becomes worse with large displacement of the SV. For high-mass LLPs (m X > 500 GeV), pair production cross sections larger than 0.07 fb are excluded for cτ 0 between 2 and 250 mm. The lowest pair production cross section excluded is 0.04 fb, at cτ 0 = 30 mm and m X > 1000 GeV. tion are calculated assuming the gluon-gluon fusion production cross section of a 125 GeV Higgs boson at 13 TeV [60]. When the long-lived scalar particle decays to a light-flavor quarkantiquark pair, branching fractions larger than 1% are excluded for cτ 0 between 1 and 340 mm with m S ≥ 40 GeV. When the long-lived scalar particle decays to a bottom quark-antiquark pair, branching fractions larger than 10% are excluded for cτ 0 between 1 and 530 mm with m S ≥ 40 GeV. These are the most stringent limits to date on this model for cτ 0 between 1 and 1000 mm. For m S = 15 GeV, where the track multiplicity of the SV is small, and the tracks are collimated due to the boost of S, the limits become worse. The limits are also worse for the case where the scalar particle decays to a bottom quark-antiquark pair, because the decays of b hadrons can produce tertiary vertices, which can be missed by the SV reconstruction we deploy in this search.
The expected and observed upper limits on the pair production cross section of long-lived gluinos in the GMSB g → g G model are shown in Fig. 9, where a branching fraction of 100% for the gluino to decay into a gluon and a gravitino is assumed. Since we do not require the reconstructed SV to have associated tracks from both jets, the two separate displaced single jets produced by the decays of the two long-lived gluinos in the g → g G model can be paired together and pass the selections, therefore the search is sensitive to the models with similar signatures. When the gluino mass is 2400 GeV, signal efficiencies are around 21, 53, and 41% in the 2017 and 2018 analysis for cτ 0 = 3, 30, and 300 mm, respectively. With the data samples collected in 2016-2018, gluino pair production cross sections larger than 0.1 fb are excluded for cτ 0 between 7 and 600 mm at m g = 2400 GeV. We then compute the upper limits on the gluino mass for different cτ 0 according to the upper limits on the pair production cross section, and a calculation at next-to-next-to-leading logarithmic precision matched to the approximated nextto-next-to-leading order predictions (NNLO approx +NNLL) of the gluino pair production cross section at √ s = 13 TeV [98][99][100][101][102][103]. Gluino masses up to 2450 GeV are excluded for cτ 0 between 6 and 550 mm. The largest gluino mass excluded is 2560 GeV with a cτ 0 of 30 mm. These limits are the most restrictive to date on this model for cτ 0 between 1 and 1000 mm. Figure 10 shows the expected and observed upper limits on the pair production cross section of the long-lived gluinos in the minisplit g → qq χ 0 1 model, assuming a branching fraction of 100% for the gluino to decay into a quark-antiquark pair and the lightest neutralino. The neutralino mass is assumed to be 100 GeV. When the gluino mass is 2400 GeV, signal efficiencies are around 31, 69, and 51% in the 2017 and 2018 analysis for cτ 0 = 3, 30, and 300 mm, respectively. With the data samples collected in 2016-2018, gluino pair production cross sections larger than 0.1 fb are excluded for proper decay lengths between 3 and 900 mm. The upper limits on the pair production cross sections are then translated into upper limits on the gluino mass for different cτ 0 , based on the NNLO approx +NNLL gluino pair production cross sections. Gluino masses up to 2500 GeV are excluded for cτ 0 between 7 and 360 mm. The largest gluino mass excluded is 2610 GeV with a cτ 0 of 30 mm. These bounds are the most stringent to date on this model for cτ 0 between 10 and 1000 mm.
The expected and observed upper limits on the pair production cross section of the long-lived gluinos in the g → tbs model are shown in Fig. 11 , where a branching fraction of 100% for the gluino to decay into top, bottom, and strange quarks is assumed. When the gluino mass is 2400 GeV, signal efficiencies are around 41, 81, and 66% in the 2017 and 2018 analysis for cτ 0 = 3, 30, and 300 mm, respectively. With the data samples collected in 2016-2018, gluino pair production cross sections larger than 0.1 fb are excluded for cτ 0 between 3 and 1490 mm at m g = 2400 GeV. We then compute the upper limits on the gluino mass for different cτ 0 according to the upper limits on the pair production cross section and the calculation of the NNLO approx +NNLL gluino pair production cross sections. Gluino masses up to 2500 GeV are excluded for cτ 0 between 3 and 1000 mm. The largest gluino mass excluded is 2640 GeV with a cτ 0 of 30 mm. These limits are the most stringent to date on this model for cτ 0 between 30 and 132 fb CMS Figure 9: Left: the 95% CL upper limits on the pair production cross section for the long-lived gluinos with m g = 2400 and 1600 GeV, where a 100% branching fraction for g → g G decays is assumed. The NNLO approx +NNLL gluino pair production cross sections for m g = 2400 and 1600 GeV, as well as their variations due to the theoretical uncertainties, are shown as horizontal lines. The solid (dashed) curves show the observed (median expected) limits, and the shaded bands indicate the regions containing 68% of the distributions of the limits expected under the background-only hypothesis. The observed limits from the CMS search utilizing the timing capabilities of the ECAL system [48] are also shown for comparison. Right: the 95% CL upper limits on the pair production cross section for the g → g G model as a function of the mean proper decay length cτ 0 and the gluino mass m g . The thick solid black (dashed red) curve shows the observed (median expected) 95% CL limit on the gluino mass as a function of cτ 0 , assuming the NNLO approx +NNLL cross sections. The thin dashed red curves indicate the region containing 68% of the distribution of the limits expected under the background-only hypothesis. The thin solid black curves represent the change in the observed limit when the signal cross sections are varied according to their theoretical uncertainties.  Figure 10: Left: the 95% CL upper limits on the pair production cross section for the long-lived gluinos with m g = 2400 GeV and 1600 GeV, where a 100% branching fraction for g → qq χ 0 1 decays is assumed. The NNLO approx +NNLL gluino pair production cross sections for m g = 2400 and 1600 GeV, as well as their variations due to the theoretical uncertainties, are shown as horizontal lines. The solid (dashed) curves show the observed (median expected) limits, and the shaded bands indicate the regions containing 68% of the distributions of the limits expected under the background-only hypothesis. Right: the 95% CL limits on the pair production cross section for the g → qq χ 0 1 model as a function of the mean proper decay length cτ 0 and the gluino mass m g . The thick solid black (dashed red) curve shows the observed (median expected) 95% CL limits on the gluino mass as a function of cτ 0 , assuming the NNLO approx +NNLL cross sections. The thin dashed red curves indicate the region containing 68% of the distribution of the limits expected under the background-only hypothesis. The thin solid black curves represent the change in the observed limit when the signal cross sections are varied according to their theoretical uncertainties.  Figure 11: Left: the 95% CL upper limits on the pair production cross section for the long-lived gluinos with m g = 2000 GeV and 1400 GeV, where a 100% branching fraction for g → tbs decays is assumed. The NNLO approx +NNLL gluino pair production cross sections for m g = 2400 and 1600 GeV, as well as their variations due to the theoretical uncertainties, are shown as horizontal lines. The solid (dashed) curves show the observed (median expected) limits, and the shaded bands indicate the regions containing 68% of the distributions of the limits expected under the background-only hypothesis. Right: the 95% CL limits on the pair production cross section for the g → tbs model as a function of the mean proper decay length cτ 0 and the gluino mass m g . The thick solid black (dashed red) curve shows the observed (median expected) 95% CL limits on the gluino mass as a function of cτ 0 , assuming the NNLO approx +NNLL cross sections. The thin dashed red curves indicate the region containing 68% of the distribution of the limits expected under the background-only hypothesis. The thin solid black curves represent the change in the observed limit when the signal cross sections are varied according to their theoretical uncertainties.
The expected and observed upper limits on the pair production cross section of the long-lived top squarks in the RPV t → b model are shown in Fig. 12, where a branching fraction of 100% for the top squark to decay into a bottom quark and a charged lepton is assumed, with equal branching fractions for e, µ, and τ. When the top squark mass is 1600 GeV, signal efficiencies are around 22, 43, and 26% in the 2017 and 2018 analysis for cτ 0 = 3, 30, and 300 mm, respectively. With the data samples collected in 2016-2018, top squark pair production cross sections larger than 0.1 fb are excluded for cτ 0 between 8 and 160 mm at m t = 1600 GeV. We then compute the upper limits on the top squark mass for different cτ 0 according to the upper limits on the pair production cross section and the calculation of the NNLO approx +NNLL top squark pair production cross sections. Top squark masses up to 1600 GeV are excluded for cτ 0 between 5 and 240 mm. The largest top squark mass excluded is 1720 GeV with a cτ 0 of 30 mm. These limits are the most stringent to date on this model in the tested cτ 0 range.
The expected and observed upper limits on the pair production cross section of the long-lived top squarks in the RPV t → d model are shown in Fig. 13, where a branching fraction of 100% for the top squark to decay into a down quark and a charged lepton, with equal branching fractions for e, µ, and τ. When the top squark mass is 1600 GeV, signal efficiencies are around 25, 48, and 29% in the 2017 and 2018 analysis for cτ 0 = 3, 30, and 300 mm, respectively. With the data samples collected in 2016-2018, and top squark pair production cross sections larger than  Figure 12: Left: the 95% CL upper limits on the pair production cross section for the long-lived top squarks with m t = 1600 GeV and 800 GeV, where a 100% branching fraction for t → b decays is assumed, with equal branching fractions for e, µ, and τ. The NNLO approx +NNLL top squark pair production cross sections for m t = 1600 and 1000 GeV, as well as their variations due to the theoretical uncertainties, are shown as horizontal lines. The solid (dashed) curves show the observed (median expected) limits, and the shaded bands indicate the regions containing 68% of the distributions of the limits expected under the background-only hypothesis. Right: the 95% CL limits on the pair production cross section for the t → b model as a function of the mean proper decay length cτ 0 and the top squark mass m t . The thick solid black (dashed red) curve shows the observed (median expected) 95% CL limits on the top squark mass as a function of cτ 0 , assuming the NNLO approx +NNLL cross sections. The thin dashed red curves indicate the region containing 68% of the distribution of the limits expected under the background-only hypothesis. The thin solid black curves represent the change in the observed limit when the signal cross sections are varied according to their theoretical uncertainties. 0.1 fb are excluded for cτ 0 between 7 and 220 mm. We then compute the upper limits on the top squark mass for different cτ 0 according to the upper limits on the pair production cross section and the calculation of the NNLO approx +NNLL top squark pair production cross sections. Top squark masses up to 1600 GeV are excluded for cτ 0 between 3 and 360 mm. The largest top squark mass excluded is 1740 GeV with a cτ 0 of 30 mm. These limits are the most restrictive to date on this model in the tested cτ 0 range.  Figure 13: Left: the 95% CL upper limits on the pair production cross section for the long-lived top squarks with m t = 1600 GeV and 800 GeV, where a 100% branching fraction for t → d decays is assumed, with equal branching fractions for e, µ, and τ. The NNLO approx +NNLL top squark pair production cross sections for m t = 1600 and 1000 GeV, as well as their variations due to the theoretical uncertainties, are shown as horizontal lines. The solid (dashed) curves show the observed (median expected) limits, and the shaded bands indicate the regions containing 68% of the distributions of the limits expected under the background-only hypothesis. Right: the 95% CL limits on the pair production cross section for the t → d model as a function of the mean proper decay length cτ 0 and the top squark mass m t . The thick solid black (dashed red) curve shows the observed (median expected) 95% CL limits on the top squark mass as a function of cτ 0 , assuming the NNLO approx +NNLL cross sections. The thin dashed red curves indicate the region containing 68% of the distribution of the limits expected under the background-only hypothesis. The thin solid black curves represent the change in the observed limit when the signal cross sections are varied according to their theoretical uncertainties.
The expected and observed upper limits on the pair production cross section of the long-lived top squarks in the dRPV t → dd model are shown in Fig. 14, where a branching fraction of 100% for the top squark to decay into two down antiquarks is assumed. When the top squark mass is 1600 GeV, signal efficiencies are around 43, 76, and 53% in the 2017 and 2018 analysis for cτ 0 = 3, 30, and 300 mm, respectively. With the data samples collected in 2016-2018, top squark pair production cross sections larger than 0.1 fb are excluded for cτ 0 between 3 and 820 mm at m t = 1600 GeV. We then compute the upper limits on the top squark mass for different cτ 0 according to the upper limits on the pair production cross section and the calculation of the NNLO approx +NNLL top squark pair production cross sections. Top squark masses up to 1600 GeV are excluded for cτ 0 between 2 and 1320 mm. The largest top squark mass excluded is 1820 GeV with a cτ 0 of 30 mm. These bounds are the most stringent to date on this model for cτ 0 between 10 and 10 4 mm.  Figure 14: Left: the 95% CL upper limits on the pair production cross section for the longlived top squarks with m t = 1600 GeV and 800 GeV, where a 100% branching fraction for t → dd decays is assumed. The NNLO approx +NNLL top squark pair production cross sections for m t = 1600 and 1000 GeV, as well as their variations due to the theoretical uncertainties, are shown as horizontal lines. The solid (dashed) curves show the observed (median expected) limits, and the shaded bands indicate the regions containing 68% of the distributions of the limits expected under the background-only hypothesis. Right: the 95% CL limits on the pair production cross section for the t → dd model as a function of the mean proper decay length cτ 0 and the top squark mass m t . The thick solid black (dashed red) curve shows the observed (median expected) 95% CL limits on the top squark mass as a function of cτ 0 , assuming the NNLO approx +NNLL cross sections. The thin dashed red curves indicate the region containing 68% of the distribution of the limits expected under the background-only hypothesis. The thin solid black curves represent the change in the observed limit when the signal cross sections are varied according to their theoretical uncertainties.

Summary
A search has been presented for long-lived particles decaying to jets, using proton-proton collision data collected with the CMS experiment at a center-of-mass energy of 13 TeV in 2017 and 2018. The results are combined with those of a previous CMS search for displaced jets using proton-proton collision data from 2016, accumulating to a total integrated luminosity of 132 fb −1 . After all selections, one event is observed in the data collected in 2017 and 2018, which is consistent with the predicted background yield. The search is designed to be model independent, and is sensitive to a large number of models predicting displaced-jets signatures with different final-state topologies.
The best current limits are set on a variety of models that have long-lived particles with mean proper decay lengths between 1 mm and 10 m. All limits are computed at the 95% confidence level. For a simplified model where pair-produced long-lived neutral particles decay to quarkantiquark pairs, pair production cross sections larger than 0.07 fb are excluded for mean proper decay lengths between 2 and 250 mm at high mass (m X > 500 GeV). For a model where the standard model-like Higgs boson decays to two long-lived scalar particles and each long-lived scalar particle decays to a down (bottom) quark-antiquark pair, branching fractions for the exotic Higgs boson decay larger than 1% (10%) are excluded for mean proper decay lengths between 1 and 340 mm (530 mm) when the scalar particle mass is larger than 40 GeV. For a supersymmetric (SUSY) model in the general gauge mediation scenario, where the long-lived gluino decays to a gluon and a lightest SUSY particle, gluino masses up to 2450 GeV are excluded for mean proper decay lengths between 6 and 550 mm. For another SUSY model in the mini-split scenario, where the long-lived gluino can decay to a quark-antiquark pair and the lightest neutralino, gluino masses up to 2500 GeV are excluded for mean proper decay lengths between 7 and 360 mm. An R-parity violating (RPV) SUSY model is also tested, where the long-lived gluino can decay to top, bottom, and strange antiquarks, and gluino masses up to 2500 GeV are excluded for mean proper decay lengths between 3 and 1000 mm. Another RPV SUSY model is studied, where the long-lived top squark can decay to a bottom quark and a charged lepton, and top squark masses up to 1600 GeV are excluded for mean proper decay lengths between 5 and 240 mm. For an RPV SUSY model, where the long-lived top squark can decay to a down quark and a charged lepton, top squark masses up to 1600 GeV are excluded for mean proper decay lengths between 3 and 360 mm. Finally, for a dynamical-RPV SUSY model, where the long-lived top squark can decay to two down antiquarks, top squark masses up to 1600 GeV are excluded for mean proper decay lengths between 2 and 1320 mm. These are the most stringent limits to date on these models.