Search for top squarks in final states with two top quarks and several light-flavor jets in proton-proton collisions at $\sqrt{s} =$ 13 TeV

Many new physics models, including versions of supersymmetry characterized by $R$-parity violation (RPV), compressed mass spectra, long decay chains, or additional hidden sectors, predict the production of events with top quarks, low missing transverse momentum, and many additional quarks or gluons. The results of a search for new physics in events with two top quarks and additional jets are reported. The search is performed using events with at least seven jets and exactly one electron or muon. No requirement on missing transverse momentum is imposed. The study is based on a sample of proton-proton collisions at $\sqrt{s} =$ 13 TeV corresponding to 137 fb$^{-1}$ of integrated luminosity collected with the CMS detector at the LHC in 2016-2018. The data are used to determine best fit values and upper limits on the cross section for pair production of top squarks in scenarios of RPV and stealth supersymmetry. Top squark masses up to 670 (870) GeV are excluded at 95% confidence level for the RPV (stealth) scenario, and the maximum observed local signal significance is 2.8 standard deviations for the RPV scenario with top squark mass of 400 GeV.


Introduction
Supersymmetry [1,2] (SUSY) is an extension of the standard model (SM) that may provide a solution to the gauge hierarchy problem [3]. In the SUSY framework, quadratically divergent radiative corrections to the Higgs boson mass parameter, dominated by loops involving the top quark, are canceled by loops with bosonic top quark superpartners (top squark, t). To avoid fine tuning, the lightest t and the superpartners of the Higgs bosons (higgsinos) must have masses near the weak scale [3][4][5][6][7][8], and could therefore have nonnegligible production cross sections at the CERN Large Hadron Collider (LHC).
Most searches for the t look for an excess of events with large missing transverse momentum p miss T originating from the undetected lightest SUSY particle (LSP) produced in t decays. It is typical in these searches to assume that the LSP is the lightest neutralino χ 0 1 , which is stable if R-parity [9] is conserved. However, it has been shown [10][11][12] that this search strategy is not sensitive to well-motivated SUSY models that predict signatures with low p miss T in models with gauge mediated SUSY breaking [13], compressed mass spectra [14,15], hidden valleys [16], or other mechanisms. As searches performed at the LHC using events with high p miss T set ever more stringent lower bounds on the t mass [17][18][19][20][21][22], searches for low-p miss T alternatives become increasingly important.
Models of R-parity violating (RPV) SUSY produce low-p miss T signatures by providing a mechanism for the LSP, in this case χ 0 1 , to decay. Among other couplings, RPV SUSY includes a trilinear Yukawa coupling between quarks and squarks that allows the χ 0 1 to decay into three quarks via an off-shell squark [9]. These couplings are typically referred to as λ ijk where i, j, and k specify the generations of the participating (s)quarks. The benchmark RPV model used in this analysis is illustrated in Fig. 1. The t decays in the typical way into a top quark and a χ 0 1 , and the χ 0 1 undergoes an RPV decay via nonzero λ 112 into three light-flavor quarks, χ 0 1 → uds. However, since this analysis does not distinguish between jets originating from quarks of the first and second generation, our results are more broadly applicable to any RPV model with coupling λ abc with a, b, c ∈ {1, 2}.
Stealth SUSY models [12,23,24] introduce a new hidden "stealth" sector of light particles with small or absent couplings to the SUSY breaking sector and finite couplings to the visible sector. Because of the weak connection to the SUSY breaking sector, SUSY is approximately conserved in the stealth sector, resulting in stealth particles that are nearly mass-degenerate with their superpartners. Production and decay of stealth particles via interactions with visible particles can be achieved through a variety of "portals" including mediation by the Higgs boson or new particles at a higher mass scale. The benchmark stealth SUSY model used in the interpretation of the results of this search (stealth SYY) [24] assumes a minimal stealth sector containing only one scalar particle S with even R-parity and its superpartner S, both of which are singlets under all SM interactions, and a portal mediated by loop interactions involving a new vectorlike messenger field (Y), the gluon (g), χ 0 1 , S, and S. Decays of the t in the stealth SYY model are illustrated in Fig. 1. Each t decays to a gluon, top quark, and S, with subsequent decays of S to S and a gravitino G and S to jets via S → gg. Because of the small mass splitting between the S and S, as well as the small G mass, the undetected G carries away very little momentum. Thus, the stealth SYY model shares the general feature of all stealth SUSY models in that it naturally produces a low-p miss T signature without R-parity violation or a special tuning of sparticle masses.
The RPV and stealth SYY models are characterized by the masses of the particles and branching fractions in the decay chain. In the benchmark RPV model, we take the χ 0 1 mass to be 100 GeV. For the benchmark stealth SYY model, the critical small S-S mass splitting is held constant at 10 GeV, and we assume a S mass of 100 GeV and a G mass of 1 GeV. For both models, a range of t masses (m t ) are considered from 300 to 1400 GeV, and all decays described above are assumed to be prompt with unity branching fractions.
In this paper, we describe a search for t pair production followed by the decay of each t into a top quark and three light-flavor jets via the benchmark RPV and stealth SYY models described above. This is the first search of its kind at the LHC. Measurements of the tt differential production cross section have been reinterpreted in the context of RPV and stealth SUSY [24, 33] and were found to yield weak constraints for the models considered in this paper.
Before describing each step in more detail in subsequent sections, we provide an overview of the analysis strategy here. The main distinguishing feature of the signals in this analysis, in addition to the presence of two top quarks, is high jet multiplicity (N jets ). The SM backgrounds arise through processes including top quark pair production (tt), multijet production from quantum chromodynamics (QCD), production of tt in association with SM weak gauge bosons or additional top quarks (tt+X), production of weak gauge bosons, and single top quark production (other). These SM processes all include additional jets from initial-and final-state radiation (ISR and FSR). The QCD background is primarily suppressed by requiring the presence of exactly one charged lepton (e or µ) arising from the leptonic decay of a top quark. Backgrounds that do not produce any top quarks are suppressed by requiring the presence of at least one jet identified as arising from the fragmentation of a bottom quark (b-tagged jet), and additionally that the invariant mass of the lepton and a b-tagged jet be consistent with the presence of a top quark.
The signal is distinguished from the dominant and irreducible tt background by means of a neural network (NN) trained to recognize differences in the spatial distribution of jets and decay kinematic distributions between signal and tt background events. Events are divided into 24 categories based on their NN score (S NN ) and N jets ; categories with higher (lower) S NN and N jets tend to be signal enriched (depleted). We perform a simultaneous fit to the number of events in data in S NN and N jets categories to estimate the total numbers of tt and potential signal events present in the data, as well as the distribution of tt events in S NN and N jets categories.
The NN output is designed to have no dependence on N jets , so that the N jets distribution of tt events can be constrained in the fit to be the same for all S NN categories. This requirement for tt N jets shape invariance is important for the analysis and will be discussed throughout the paper. This paper is organized as follows. We introduce the CMS detector and methods for event reconstruction and selection in Section 2. Samples of simulated events are described in Section 3. The estimation and modeling of SM backgrounds are explained in Section 4, and the description of the treatment of systematic uncertainties is in Section 5. Finally, the results and their interpretation are in Section 6, followed by the summary in Section 7.

Experimental techniques
The search is performed using a data sample of proton-proton (pp) collisions at √ s = 13 TeV, corresponding to an integrated luminosity of 137 fb −1 , collected in 2016-2018 with the CMS detector at the LHC. Data and simulation samples from four periods (2016, 2017, 2018A, 2018B) are treated separately in order to address variations in detector and LHC conditions. Data from 2018 are divided into two samples (2018A and 2018B), with 2018B corresponding to the period when a detector malfunction prevented readout from 3% of the hadron calorimeter. In this section, we define reconstructed physics objects and describe the selection criteria for events in the signal region (SR) and the control region (CR) of the analysis.
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, and a brass and scintillator hadron calorimeter, 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 kinematic variables, can be found in Ref.
[34]. The CMS trigger system is described in Ref. [35]. Events are selected using triggers that require the presence of at least one electron or one muon. The minimum transverse momentum p T threshold is 27 (35) GeV for electrons and 24 (24) GeV for muons in 2016 (2017-2018). The triggers at these thresholds require the lepton to be isolated from tracks and calorimeter deposits in the detector. Events may also be selected from single-lepton triggers with higher p T thresholds, 115 GeV for electrons and 50 GeV for muons, with no isolation requirements. The combined trigger efficiency varies from 80% for leptons with p T close to the lower thresholds to greater than 95% for leptons with p T > 120 GeV.
Events are reconstructed using the particle-flow (PF) algorithm [36], which reconstructs particles in an event using an optimized combination of information from the various elements of the CMS detector and identifies each as a photon, electron, muon, charged hadron, or neutral hadron. These particles are further clustered into jets as described below.
The reconstructed vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex, where the physics objects are the jets, clustered using the antik T algorithm [37, 38] with the charged-particle tracks assigned to the vertex as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets [39]. Charged-particle tracks associated with vertices from other pp interactions (pileup) are removed from further consideration. The primary vertex is required to lie within 24 cm of the interaction point along the beam axis, and within 2 cm in the plane transverse to the beam axis.
Electrons and muons must satisfy p T > 30 GeV and |η| < 2.4. For the analysis of the 2017 and 2018 data, the electron p T threshold is increased to 37 GeV to account for the higher trigger threshold. The lepton identification requirements are the "tight" criteria for electrons [40] and the "medium" criteria for muons [41]. Leptons must be isolated within a cone of radius R = √ (∆φ) 2 + (∆η) 2 that scales as 1/p T between a maximum of 0.2 for leptons with p T < 50 GeV and a minimum of 0.05 for lepton p T > 200 GeV [42].
Jets are clustered from the reconstructed PF particles using the anti-k T algorithm with a distance parameter of 0.4. Criteria are applied to remove events with jets arising from instrumental effects or reconstruction failures [43,44]. The reconstructed jet energies are corrected for the nonlinear response of the detector [45,46] and for contributions from neutral hadrons from pileup [47]. Jets are required to have p T > 30 GeV and |η| < 2.4. Jets overlapping with a selected lepton within a cone of radius R = 0.4 are removed. A neural network-based algorithm [48] is used to identify b quark jets; for jets with p T around 30 GeV, the algorithm has an efficiency of 65% and a misidentification rate for light-flavor jets (including gluon jets) of 1%.
In addition to the trigger and vertex criteria above, events in the SR must contain exactly one isolated electron or muon and at least seven jets, at least one of which should be b tagged. Samples with seven and eight jets include a small number of expected signal events, but are included in the SR to constrain the background. To further reject the QCD background, we require the scalar sum of jet p T (H T ) to exceed 300 GeV. To suppress non-tt backgrounds, we require the invariant mass of the system formed by the b-tagged jet and the lepton to be between 50 and 250 GeV. If there is more than one b-tagged jet in the event, the invariant mass of each b-tagged jet and the lepton is considered, and at least one combination is required to meet the above criterion. No requirement is made on the event p miss T . In addition to the SR, a signal-depleted control region (CR) dominated by QCD background is defined with the dual purpose of determining the QCD contribution to the SR and verifying the important assumption of tt N jets shape invariance with S NN . Despite being dominated by QCD background, the CR is useful for confirming tt N jets shape invariance because many of the jets used as inputs to the NN arise from QCD radiation, which is common to the tt and QCD backgrounds; this claim is verified in Section 5. The CR is defined similarly to the SR with the differences being that the lepton is required to be a muon; the muon is required to fail the SR isolation requirement; there is no requirement for a b-tagged jet, nor on the invariant mass of the lepton and b-tagged jet; the only trigger used is the high-threshold muon trigger without an isolation requirement; and the muon p T threshold is 55 GeV.

Simulated event samples
Simulated event samples are used in the estimation of the expected number of SM background and signal events passing the SR selection. Top quark pair and single top quark events produced in the t channel are generated with the next-to-leading-order (NLO) POWHEG v2.0 [49][50][51][52][53] generator, while single top quark events in the tW channel are generated with POWHEG v1.0 [52]. Single top quark production in the s channel, as well as rare SM processes such as ttZ and ttW are generated at NLO accuracy with the MADGRAPH5 aMC@NLO v2.2.2 program. The MADGRAPH5 aMC@NLO v2.2.2 generator [54,55] is used in the leading-order (LO) mode to simulate QCD and W+jets events.
For the signal, top squark pair production events are generated using MADGRAPH5 aMC@NLO in LO mode, including up to two additional partons in the matrix element calculation. The top squarks are decayed using PYTHIA v8.212 (2016) or 8.226 (2017-2018) [56] according to the signal models described in Section 1. The signal production cross section (σ t t ) is calculated as a function of m t using approximate next-to-NLO (NNLO) plus next-to-next-to-leading-logarithm (NNLL) calculations [57,58].
The generation of these processes is based on either LO or NLO parton distribution functions (PDFs) using NNPDF3.0 [59] for the simulated samples corresponding to 2016 detector conditions, and using the NNLO PDF sets from NNPDF3.1 [60] for the 2017 and 2018 simulated samples. Parton showering and hadronization are simulated with PYTHIA using underlying event tune CUETP8M1 [61] for 2016 samples, except for tt production which used tune CUETP8M2T4 [62], or PYTHIA with tune CP5 (CP2) [63] for all 2017 and 2018 background (signal) samples. To model the effects of pileup, simulated events are generated with a nominal distribution of pp interactions per bunch crossing and then reweighted to match the corresponding distribution in data. The CMS detector response is simulated using a GEANT4-based model [64], and event reconstruction is performed in the same manner as for collision data. The most precise cross section calculations available are used to normalize the SM simulated samples, corresponding to NLO or NNLO accuracy in most cases [54,[65][66][67][68][69][70][71].
The simulation is corrected to eliminate small discrepancies between data and simulation in the trigger efficiency, lepton selection efficiency, and b tagging efficiency. Analysis-specific corrections for the H T distribution in tt simulation, parameterized as functions of N jets and H T , are obtained in a signal-depleted sample identical to the SR, except for the requirement 5 ≤ N jets ≤ 7. Events with N jets = 7 are common to the SR, but as mentioned above, this sample has low signal contamination. The correction is parameterized with an exponential function in H T with parameters depending linearly on N jets in order to extend the correction into the N jets > 7 SR. The H T correction is small at low H T and 20-40% at H T = 1500 GeV, depending on N jets .

Background estimation
Simulated background events passing the SR selection requirements predominantly originate from tt production, with contributions of less than 10% from QCD, and a few percent from the remaining minor backgrounds including tt production in association with a vector boson, single top quark production, and W+jets.
As introduced in Section 1, the crux of the analysis is to estimate the dominant tt background in four bins of S NN and six N jets bins using a simultaneous binned maximum-likelihood fit constraining the tt N jets shape to be the same in all S NN categories. Event yields, as well as the N jets and S NN distributions, are fixed at values determined from a signal-depleted data control sample for the QCD background and from simulation for the minor backgrounds, as described later in this section. The yield and N jets shape of the tt background, along with the signal strength, are determined in the fit; signal strength is defined as the ratio of the fit signal event yield to the one predicted by SUSY.
The NN is trained to discriminate between signal and tt background by exploiting differences in the event shape and distributions of the kinematic variables. The gradient reversal technique [72] is used to minimize dependence of the NN output on N jets , as required by the primary assumption that the tt N jets shape is the same in all S NN categories. All NN input variables are computed in an approximate center-of-mass frame defined by all jets in the event with p T > 30 GeV and |η| < 5. The NN input variables include the four-vector components for the seven jets in the event with the highest momentum in the center-of-mass frame, the four-vector components of the lepton in the event, the second through fifth Fox-Wolfram moments [73] normalized by the first moment, and the three eigenvalues of the sphericity tensor [74] normalized by the sum of the eigenvalues. The Fox-Wolfram moments and sphericity tensor eigenvalues, which are computed from the same seven highest momentum jets, quantify the distribution of jet energy in the event, which tends to be more spherical for signal t pair production than for the tt background.
For the NN training, simulated tt events are used for the background sample, and a mixture of RPV and stealth SYY simulated events with m t from 350-850 GeV is used as the signal sample. In this way, the NN can identify common features among all signal samples ensuring a search with broad sensitivity. Reflecting differences in simulation between the data taking periods, as described in Section 3, a single training is used for 2017, 2018A, and 2018B, with a separate training used for 2016. The S NN distributions for the simulated background, several signal models, and the 2016 and 2017+2018 data are shown in Fig. 2. Events Events For each of the six N jets bins, events are divided into four S NN bins: S NN,1 (lowest S NN ), . . . , S NN,4 (highest S NN ). The S NN bin boundaries are chosen separately for each N jets bin such that the expected significance for the 550 GeV RPV signal model, which has expected significance close to 5 standard deviations (σ), is maximized, under the constraint that the fraction of simulated tt events in each S NN bin is the same for all N jets bins. For example, the fraction of all events in each N jets bin falling into the S NN,1 bin is constrained to be approximately 56%, while the fraction of events falling into the S NN,4 bin is constrained to be approximately 2.4%. This constraint removes the small dependence of S NN on N jets that remains after NN training with gradient reversal.
In the maximum-likelihood fit, the tt N jets distribution is parameterized with a function inspired by QCD jet scaling patterns [75] in which the ratio R(i) = M i+1 /M i , where M i is the number of events with N jets = i, can be described by a falling "Poisson" component at low N jets and a constant "staircase" component at high N jets . This ratio is well modeled by the function Notice that a 0 = f (7), a 1 = f (9), and a 2 is the asymptotic value for large i. This particular parameterization was chosen to avoid large correlations between the fit parameters. The N jets distribution for each S NN bin j (see Fig. 4) is modeled using a recursive expression given by are normalization parameters that are floating in the fit. The last N jets bin considered is an inclusive N jets ≥ 12 bin, such that i ∈ [7,12]. In the maximumlikelihood fit, the free parameters consist of the three shape parameters a 0 , a 1 , and a 2 ; the four normalization parameters Y j 7 ; the signal strength; and all nuisance parameters related to systematic uncertainties described in Section 5.
The QCD background yield parameters are fixed in the fit at the values determined from the CR. More specifically, the QCD background prediction for each N jets -S NN bin in the SR is given by the yield for the same bin in the CR in data, after subtraction of the non-QCD backgrounds as predicted from simulation, multiplied by the ratio of SR to CR yields in simulation (R QCD ). This procedure is verified with a closure test in the simulation. The yield parameters from the minor backgrounds are also kept fixed in the fit at the values predicted by simulation. While the yield parameters are fixed in the fit, the ultimate contributions from QCD and minor backgrounds vary according to the constrained nuisance parameters related to systematic uncertainties in those fit components.

Systematic uncertainties and fit validation
As described in Section 4, an unbiased estimate of the dominant tt background is obtained from the fit to data as long as the tt N jets shape is the same for all four S NN bins. By construction, N jets shape invariance is achieved in the simulation with an N jets -specific S NN binning as described in the previous section. Thus, systematic uncertainties on the tt background are important to the degree that they violate the assumption that the S NN binning determined in simulation also applies to the data. We quantify how each source of uncertainty causes deviations from the assumed N jets shape invariance by comparing the nominal N jets shape to the N jets shapes in all S NN bins after performing the relevant systematic variation to the tt simulation. Each systematic variation is associated with a constrained nuisance parameter in the fit. The deviation in shape for each N jets distribution, derived from the ratio of the post-variation shape divided by the nominal shape, changes linearly with the associated nuisance parameter for the systematic variation, while preserving the normalization of the distribution.
Sources of tt shape uncertainty include uncertainty in aspects of event generation including PDFs, choice of renormalization and factorization scales (µ R , µ F scales), and parton shower modeling, which is itself composed of aspects related to modeling of ISR, FSR, color reconnection in the parton shower, matrix element-parton shower matching scale (ME-PS), underlying event (UE tune), and pileup modeling. The uncertainty due to the choice in (µ R , µ F ) scales is determined by independently varying both by factors of 2.0 and 0.5 excluding the variations (2.0, 0.5) and (0.5, 2.0) [55,76,77]. The ISR and FSR uncertainties originate from variations of the renormalization scale for the parton shower by factors 0.5 and 2.0, effectively varying the value of α S . The color reconnection uncertainty is calculated by allowing resonant decays to occur before the merging of multi-parton systems. The ME-PS uncertainty is obtained by varying the POWHEG parameter that governs ME-PS matching about its nominal value according to h damp = 1.379 +0.926 −0.505 times the top quark mass [63]. The UE tune uncertainty comes from variation of the PYTHIA parameters that control the modeling of the underlying event as described in Ref. [63]. The total inelastic pp cross section is changed by 5% to estimate the uncertainty related to pileup [78].
Sources of tt shape uncertainty related mostly to aspects of detector simulation include determination of jet energy scale (JES) and resolution (JER), modeling of the b tagging efficiency, modeling of the efficiency for lepton triggers, identification, and isolation (lepton efficiencies); residual mismodeling of H T , jet p T , and jet mass; and use of the CR for measuring deviations from the assumption of N jets shape invariance.
The uncertainty in the modeling of H T in the tt simulation is composed of four separate components. The first H T uncertainty (primary) is taken as the full difference in the tt background shape with and without the H T correction. The second H T uncertainty (validation) is taken as the difference between the simulation with nominal H T correction (described in Section 3) and the observed H T distribution in the signal-depleted SR sample with N jets = 8. The third and fourth H T uncertainties address the choices of parameterization of the H T correction as functions of H T and N jets . For these, we take the uncertainty as the difference between the nominal correction and two alternate corrections that use the H T = 2000 GeV correction for all events with H T > 2000 GeV (H T -parameterization) and the N jets = 7 correction for all values of N jets (N jets -parameterization).
Comparisons of data and simulation in the CR show that the simulation predicts distributions with higher values of jet p T and mass than observed. The observed discrepancy at jet p T (mass) of 400 (50) GeV depends on jet p T rank and is small for the highest p T jet in each event growing to approximately 50% for the jet with sixth-highest p T in each event. Similar trends are observed in the signal-depleted, tt-dominated SR with N jets = 7. In the CR, the discrepancy in the falling tail of each distribution is minimized when the p T (mass) of each jet is scaled by the value 0.95, 0.95, 0.95, 0.95 (0.95, 1.01, 0.98, 0.98) for 2016, 2017, 2018A, and 2018B, respectively. Thus, the related tt shape uncertainty is taken to be the resulting difference between scaled and nominal simulated tt distributions. The dependence on jet p T rank indicates that the discrepancy arises predominantly in the event generation; however, we choose to estimate the associated systematic uncertainty separately for each data taking period to include potential effects of detector response simulation. The H T correction is omitted from the determination of these jet p T and mass uncertainties to avoid double counting of H T mismodeling effects. In addition, because the estimation of jet p T and mass uncertainties relies on variable scaling (rather than event reweighting), the uncertainties include effects of changes in the S NN for each event, which is not included in the H T uncertainty.
As mentioned above, the use of N jets -dependent S NN binning ensures that the N jets shape is the same in all four S NN bins in simulation, and the use of the same binning in the data assumes that the N jets -S NN dependence is well modeled in the simulation. This assumption is confirmed and a related systematic uncertainty is determined by comparing the N jets shapes (in five uniform S NN bins) for data and simulation in the CR. , and the data in the QCD background-dominated CR (Data CR ). Agreement between QCD MC CR and tt MC SR demonstrates that QCD background-dominated data in the CR are a good proxy for tt-dominated data in the SR, and agreement between QCD MC CR and Data CR verifies that the dependence of the N jets shape on S NN is well modeled in the simulation. Similar agreement is found for the R M distributions for the other N jets bins and data periods. The uncertainty related to the combination of both effects is taken as the difference between tt MC SR and Data CR . For the QCD background, the shape is obtained from data in the CR, and the normalization is set with R QCD . Because the systematic uncertainties in the simulation largely cancel in the R QCD ratio, the uncertainty in R QCD is dominated by the statistical uncertainty of simulated samples and ranges from 15-25% depending on data collection period.
Sources of systematic uncertainty in the predictions for signals and the minor backgrounds include PDFs, JES, JER, b tagging efficiency, lepton efficiency, trigger efficiency, (µ R , µ F ) scales, cross sections for the minor backgrounds, and a 2.3-2.5% uncertainty in the integrated luminosity [79][80][81]. Since the signal and minor backgrounds are estimated directly from simulation, related uncertainties are included as the full effect of the systematic variation on the yields in each N jets and S NN bin, thereby taking into account normalization effects as well as shape changes.
Uncertainties derived from comparisons of data and simulation separately in each data taking period (related to pileup, JES, JER, b tagging efficiency, lepton efficiencies, H T corrections, N jets shape invariance, and integrated luminosity) are treated as uncorrelated among all data samples. Uncertainties related to parton shower modeling are treated as fully correlated for 2017, 2018A, and 2018B, while the corresponding uncertainties for 2016 are uncorrelated with those from the other data taking periods; uncertainties related to (µ R , µ F ) scales and cross sections for the minor backgrounds are treated as correlated between all four periods. Table 1 shows the impact of the systematic uncertainties on the expected event yields for the tt background, minor backgrounds, and the RPV signal model with m t = 550 GeV. For sources of uncertainty for which the size of the impact depends on N jets and S NN , a representative range of values is listed along with the maximum value from all bins.

Results and interpretation
The results of the fit to 2016, 2017, 2018A, and 2018B data sets with the signal strength fixed to zero (background-only fit) are shown along with the observed number of events in Fig. 4; each column (row) in the plot array corresponds to a specific S NN bin (data set). The expected distributions for top squark pair production in the specific RPV (m t = 450 GeV) and stealth SYY models (m t = 850 GeV) described in Section 1 are overlaid for illustration purposes. The lower panel of each plot displays the difference between the observed number of events and the total number of expected events determined by the fit divided by the statistical uncertainty associated with the observed number of events (δ) as black points with error bars denoting δ. The   blue band shows the total uncertainty in the fit determined from the full fit covariance matrix in order to account for the correlations among fit parameters. Figure 5 shows the results of the same background-only fit summed over S NN bins and data periods with separate contributions from each background.
The data are also used to determine the 95% confidence level (CL) upper limits on σ t t and the signal strength p-values [82] for both the RPV and stealth SYY models obtained using the CL s approach [83][84][85] with asymptotic formulae [86] and the profile likelihood ratio as the test statistic. Figure 6 shows the expected and observed cross section limits as a function of m t for the benchmark RPV and stealth SYY signal models. Comparing to the predicted cross section, these limits correspond to the exclusion of top squark masses in the range 300-670 and 300-870 GeV for the benchmark RPV and stealth SYY models, respectively. Figure 7 shows the local p-value [82] of the signal strength, as a function of m t , obtained from fits to the data with each signal strength as a free parameter for both the RPV and stealth SYY models. The p-value quantifies the probability for the background to produce an upward fluctuation at least as large as that observed. Fits are performed and p-values obtained separately for each data set, as well as in a simultaneous fit to all data sets. We observe the most extreme p-value to be 0.003, which corresponds to a local significance of 2.8 σ and a best fit signal strength of 0.21 ± 0.07 for the RPV model with m t = 400 GeV assuming unity branching fractions for the decays described in Section 1.
The 2.8 σ local significance for the RPV model with m t = 400 GeV is understood to arise from a combination of two effects. First, although the level of agreement between the observed data and the background expectation shown in Fig. 4 is reasonable, the agreement improves when the signal is included in the fit, contributing approximately 1.1 σ to the significance. Second, the constrained nuisance parameters (NP) are pulled less from their initial values when the signal  [pb] t t σ 95% CL upper limit on   is included in the fit, contributing approximately 1.7 σ to the significance. This second effect is illustrated in Fig. 8 which shows for each of a selection of NP : the fit value (θ) and uncertainty (δ θ ) from both the background-only fit (b) and the signal+background fit (s+b), as well as the ∆χ 2 ≡ χ 2 (s + b) − χ 2 (b) difference of χ 2 ≡ (θ/δ θ ) 2 from the two fits. A θ value of one indicates that the fit value of the nuisance parameter is one standard deviation from its nominal value, and a δ θ value less than one shows that the uncertainty is reduced in the fit relative to its initial value. All NPs have θ values below one for the background-only fit, and several NPs related to tt modeling are constrained with δ θ in the range of 0.25-0.40. Figure 8 also shows the cumulative and total sums of ∆χ 2 for the NPs, with the sum for all NP of ∑ ∆χ 2 = −3.0 corresponding to an approximate contribution to the signal significance of | ∑ ∆χ 2 | = 1.7 σ.

Summary
A first of its kind search for top squark pair production with subsequent decay characterized by two top quarks, additional gluons or light-flavor quarks, and low missing transverse momentum (p miss T ) is described. Events containing exactly one electron or muon and at least seven jets, of which at least one should be b tagged, are selected from a sample of proton-proton collisions at √ s = 13 TeV corresponding to an integrated luminosity of 137 fb −1 collected with the CMS detector in 2016-2018. No requirement is made on p miss T . The dominant tt background is predicted from data using a simultaneous fit of the jet multiplicity distribution across four bins of a neural network score.
The results are interpreted in terms of top squark pair production in the context of R-parity violating (RPV) and stealth supersymmetry models. Top squark masses (m t ) up to 670 GeV are excluded at 95% confidence level for the RPV model in which the top squark decays to a top quark and the lightest neutralino, which subsequently decays to three light-flavor quarks via an off-shell squark through a trilinear coupling λ . Top squark masses up to 870 GeV are excluded for the stealth supersymmetry model in which the top squark decays to a top quark, three gluons, and a gravitino via intermediate hidden sector particles. The maximum observed local significance is 2.8 standard deviations corresponding to a best fit signal strength of 0.21 ± 0.07 for the RPV model with m t = 400 GeV.

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 centers and personnel of the Worldwide LHC Computing Grid and other centers 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, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies:   [22] ATLAS Collaboration, "Search for top-squark pair production in final states with one lepton, jets, and missing transverse momentum using 36 fb −1 of √ s = 13 TeV pp collision data with the ATLAS detector", JHEP 06 (2018) 108, doi:10.1007/JHEP06(2018)108, arXiv:1711.11520.   [40] CMS Collaboration, "Performance of electron reconstruction and selection with the CMS detector in proton-proton collisions at √ s = 8 TeV", JINST 10 (2015) P06005, doi:10.1088/1748-0221/10/06/P06005, arXiv:1502.02701.