Search for squarks and gluinos in events with an isolated lepton, jets, and missing transverse momentum at √s=13 TeV with the ATLAS detector

The results of a search for squarks and gluinos in final states with an isolated electron or muon, multiple jets and large missing transverse momentum using proton-proton collision data at a center-of-mass energy of ﬃﬃﬃ s p ¼ 13 TeV are presented. The data set used was recorded during 2015 and 2016 by the ATLAS experiment at the Large Hadron Collider and corresponds to an integrated luminosity of 36 . 1 fb − 1 . No significant excess beyond the expected background is found. Exclusion limits at 95% confidence level are set in a number of supersymmetric scenarios, reaching masses up to 2.1 TeV for gluino pair production and up to 1.25 TeV for squark pair production.


Introduction
Supersymmetry (SUSY) [1][2][3][4][5][6] is a theoretical framework of physics beyond the Standard Model (SM) which predicts for each SM particle the existence of a supersymmetric partner (sparticle) differing by half a unit of spin. The partner particles of the SM fermions (quarks and leptons) are the scalar squarks (q) and sleptons (˜ ). In the boson sector, the supersymmetric partner of the gluon is the fermionic gluino (g), whereas the supersymmetric partners of the Higgs (higgsinos) and the electroweak gauge bosons (winos and bino) mix to form charged mass eigenstates (charginos) and neutral mass eigenstates (neutralinos). In the minimal supersymmetric extension of the Standard Model (MSSM) [7,8] two scalar Higgs doublets along with their higgsino partners are necessary, resulting in four chargino states (χ ± 1,2 ) and four neutralinos (χ 0 1,2,3,4 ). SUSY addresses the SM hierarchy problem [9][10][11][12] provided that the masses of at least some of the supersymmetric particles (most notably the higgsinos, the top squarks and the gluinos) are near the TeV scale.
In R-parity-conserving SUSY [13], gluinos or squarks are pair produced at the Large Hadron Collider (LHC) via the strong interaction and decay either directly or via intermediate states to the lightest supersymmetric particle (LSP). The LSP, which is assumed to be the lightest neutralino (χ 0 1 ) in this paper, is stable and weakly interacting, making it a candidate for dark matter [14,15].
The decay topologies targeted in this paper are largely inspired by decay chains that could be realized in the pMSSM scenario, which is a two-dimensional subspace of the 19-parameter phenomenological Minimal Supersymmetric Standard Model (pMSSM) [16,17]. Four SUSY models with gluino or squark pair production and different decay topologies are considered. The first two models, referred to as the gluino and squark one-step models for the rest of this paper, are SUSY simplified models [18][19][20] in which pair-produced gluinos or squarks decay via the lightest chargino (χ ± 1 ) to the LSP. In the model with gluino production, the gluino decays to the lightest chargino and two SM quarks viag → qq χ ± 1 , as illustrated in Figure 1 (left). The gluino decay is assumed to proceed via virtual first-and second-generation squarks, hence no bottom or top quarks are produced in the simplified model. The chargino then decays to the LSP by emitting an on-or off-shell W boson,χ ± 1 → W ( * )±χ 0 1 , depending on the available phase space. In the MSSM this decay chain is realized when the gluino decays, via a virtual squark that is the partner particle of the left-handed SM quark, to the chargino with a dominant wino component. In the squark production model, the squark decays to the chargino viaq → q χ ± 1 , followed by the same chargino decay, as illustrated in Figure 1 (middle). The third model, referred to as the gluino two-step model for the rest of this paper, assumes gluino pair production with a subsequent decay to the chargino viag → qq χ ± . The chargino then decays via emission of an on-or off-shell W boson to the second lightest neutralino according toχ ± → W ±χ 0 2 . In the last step of the cascade, the second lightest neutralino decays via emission of a Z boson to the LSP. The decay chain of this signal model is illustrated in Figure 1 (right). The model is used as a proxy for SUSY scenarios with many decay products in the final state. Within the MSSM, additional decay modes lead to a significant reduction in the cross-section times branching fraction for this particular decay.
Finally, the fourth set of SUSY models, the pMSSM model, is selected to have a bino-dominated neutralino as the LSP, kinematically accessible gluinos, and a higgsino-dominated multiplet at intermediate mass.
The higgsino multiplet contains two neutralinos (theχ 0 2 andχ 0 3 ) and a chargino. The decays proceed predominantly via virtual third-generation supersymmetric quarks due to their enhanced couplings with the higgsinos. Examples of dominant characteristic decay chains of this model for mχ ± 1 ∼ < 500 GeV and mg ∼ > 1200 GeV areg → ttχ 0 2,3 andg → tbχ ± 1 , withχ 0 2,3 decaying to Z/hχ 0 1 andχ ± 1 to W ±χ0 1 . In this search, the experimental signature consists of a lepton (electron or muon), several jets, and missing transverse momentum (E miss T ) from the undetectable neutralinos and neutrino(s). Depending on the sparticle masses of the model considered, different amounts of energy are available in their decays. Therefore, the number of leptons and jets in the final state, as well as their kinematic properties, depend on the mass spectrum in the model of interest. Four signal regions with jet multiplicities ranging from two to six are defined to provide sensitivity to a broad range of mass spectra in the gluino and squark one-step models. For the two-step and pMSSM models, a dedicated signal region requiring nine jets is constructed to take advantage of the large jet multiplicities in these models. In each signal region, the event yield is compared with the SM prediction, which is estimated using a combination of simulation and observed data in control regions. The results of all Run-1 ATLAS searches targeting squark and gluino pair production are summarized in Ref. [28]. The same SUSY models considered in this paper were also targeted in other Run-2 ATLAS searches using different experimental signatures [29][30][31].
This paper is structured as follows. After a brief description of the ATLAS detector in Section 2, the simulated data samples for the background and signal processes used in the analysis as well as the dataset and the trigger strategy are detailed in Section 3. The reconstructed objects and quantities used in the analysis are described in Section 4 and the event selection is presented in Section 5. The background estimation and the systematic uncertainties associated with the expected event yields are discussed in Sections 6 and 7, respectively. Finally, the results of the analysis are presented in Section 8, and are followed by a conclusion. ATLAS [32] is a general-purpose detector with a forward-backward symmetric design that provides almost full solid angle coverage around the interaction point. 1 The main components are the inner detector (ID), which is surrounded by a superconducting solenoid providing a 2 T axial magnetic field, the calorimeter system, and the muon spectrometer (MS), which is immersed in a magnetic field generated by three large superconducting toroidal magnets. The ID provides track reconstruction within |η| < 2.5, employing pixel detectors close to the beam pipe, silicon microstrip detectors at intermediate radii, and a straw-tube tracker with particle identification capabilities based on transition radiation at radii up to 1080 mm. The innermost pixel detector layer, the insertable B-layer [33], was added during the shutdown between LHC Run 1 and Run 2, at a radius of 33 mm around a new, narrower, beam pipe. The calorimeters cover |η| < 4.9. The forward region (3.2 < |η| < 4.9) is instrumented with a liquid-argon (LAr) calorimeter for both the electromagnetic and hadronic measurements. In the central region, a lead/LAr electromagnetic calorimeter covers |η| < 3.2, while the hadronic calorimeter uses two different detector technologies, with scintillator tiles (|η| < 1.7) or liquid argon (1.5 < |η| < 3.2) as the active medium. The MS consists of three layers of precision tracking chambers providing coverage over |η| < 2.7, while dedicated fast chambers allow triggering over |η| < 2.4. The ATLAS trigger system used for real-time event selection [34] consists of a hardware-based first-level trigger and a software-based high-level trigger.

Simulated event samples and data samples
Three simplified SUSY signal models and a set of pMSSM scenarios are considered in this search. Gluinos or squarks are assumed to be produced in pairs (gg orqq). In the case of the simplified models, 100% branching ratios to the decay of interest are assumed.
The gluino/squark one-step simplified models have three free parameters: the masses of the gluino or squark (mg /q ), the lightest chargino (mχ ± 1 ), and the lightest neutralino (mχ 0 1 ). Other sparticles that do not appear in the decay chain are set to have a high mass. To probe a broad range of SUSY mass spectra, two model parameterizations are considered. In the first type, mg /q and the mass ratio x ≡ (mχ ± 1 − mχ 0 1 )/(mg /q − mχ 0 1 ) are free parameters, while mχ 0 1 is fixed to 60 GeV. In the second type, mg /q and mχ 0 1 are free parameters, while mχ ± 1 is fixed by setting x = 1/2. For the rest of this paper, the former type is referred to as variable-x and the latter one is referred to as x = 1/2.
The gluino two-step simplified model has two free parameters that are varied to probe different mass configurations: the masses of the gluino (mg) and the lightest neutralino (mχ 0 1 ). The masses of the lightest chargino and the second-lightest neutralino are constrained to be mχ ± 1 = (mg + mχ 0 1 )/2 and mχ 0 2 = (mχ ± 1 + mχ 0 1 )/2, respectively. All other sparticles are kinematically inaccessible. In the pMSSM scenario, the sparticle masses are varied by scanning the gluino mass parameter M 3 (related to mg) and the bilinear Higgs mass parameter µ (related to mχ ± 1 and mχ 0 2 ). The scan ranges are 690 GeV < M 3 < 2140 GeV and −770 GeV < µ < −160 GeV. The bino mass parameter M 1 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector and the z-axis along the beam pipe. The x-axis points from the IP to the center of the LHC ring, and the y-axis points upward. Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis. The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2). Rapidity is defined as y = 0.5 ln[(E + p z )/(E − p z )] where E denotes the energy and p z is the component of the momentum along the beam direction. (related to mχ 0 1 ) was set to 60 GeV. The remaining model parameters, defined in Ref. [35], are set to TeV, such that the mass of the lightest Higgs boson is compatible with 125 GeV and all other sparticles are kinematically inaccessible. Mass spectra consistent with electroweak symmetry breaking were generated using SOFTSUSY 3.4.0 [36] and the decay branching ratios were calculated with SDECAY/HDECAY 1.3b/3.4 [37].
The signal samples were generated at leading order (LO) using M G 2.2.2 [38] with up to two extra partons in the matrix element, interfaced to P 8.186 [39] for parton showers and hadronization. The CKKW-L matching scheme [40] was applied for the matching of the matrix element and the parton shower, with a scale parameter set to a quarter of the mass of the sparticle produced. The ATLAS A14 [41] set of tuned parameters (tune) was used for the shower and the underlying event, together with the NNPDF2.3 LO [42] parton distribution function (PDF) set. The E G 1.2.0 program [43] was used to describe the properties of the bottom and charm hadron decays in the signal samples.
The signal cross-sections were calculated at next-to-leading order (NLO) in the strong coupling constant, adding the resummation of soft gluon emission at next-to-leading-logarithmic accuracy (NLL) [44][45][46][47][48]. The nominal cross-section and its uncertainty are taken from an envelope of cross-section predictions using different PDF sets and factorization and renormalization scales, as described in Ref. [49], considering only the four light-flavor left-handed squarks (ũ L ,d L ,s L , andc L ).
The simulated event samples for the signal and SM backgrounds are summarized in Table 1. Additional samples are used to assess systematic uncertainties, as explained in Section 7.
To generate tt and single-top-quark events in the Wt and s-channel [50], the P -B v2 [51] event generator with the CT10 [52] PDF set in the matrix-element calculations was used. Electroweak t-channel single-top-quark events were generated using the P -B v1 event generator. This event generator uses the four-flavor scheme for the NLO matrix-element calculations together with the fixed four-flavor PDF set CT10f4. For all top quark processes, top quark spin correlations are preserved (for the singletop t-channel, top quarks are decayed using MadSpin [53]). The parton shower, fragmentation, and the underlying event were simulated using P 6.428 [54] with the CTEQ6L1 [55] PDF set and the corresponding P 2012 tune (P2012) [56]. The top quark mass was set to 172.5 GeV. The EvtGen 1.2.0 program was also used to describe the properties of the bottom and charm hadron decays in the tt and the single-top-quark samples. The h damp parameter, which controls the p T of the first additional emission beyond the Born configuration, was set to the mass of the top quark. The main effect of this is to regulate the high-p T emission against which the tt system recoils. The tt events are normalized using the cross-sections computed at next-to-next-to-leading order (NNLO) with next-to-next-to-leadinglogarithmic (NNLL) corrections [57]. The single top quark events are normalized using the NLO+NNLL cross-sections for the Wt-channel [58] and to the NLO cross-sections for the tand s-channels [59].
Events containing W or Z bosons with associated jets (W/Z+jets) [60] were simulated using the S 2.2.1 event generator [61]. Matrix elements were calculated for up to two partons at NLO and four partons at LO using the Comix [62] and OpenLoops [63] generators. They were merged with the S 2.2.1 parton shower [64] with massive band c-quarks using the ME+PS@NLO prescription [65]. The NNPDF3.0 NNLO PDF set [66] was used in conjunction with a dedicated parton shower tuning developed by the S authors. The W/Z+jets events are normalized using their NNLO cross-sections [67].
The diboson samples [68] were generated using the S 2.1.1 and 2.2.1 event generators using the CT10 and NNPDF3.0 PDF sets, respectively. The fully leptonic diboson processes were simulated including final states with four charged leptons, three charged leptons and one neutrino, two charged leptons and two neutrinos, and one charged lepton and three neutrinos. The semileptonic diboson processes were simulated with one of the bosons decaying hadronically and the other leptonically. The processes were calculated for up to one parton (for Z Z) or no additional partons (for WW, W Z) at NLO and up to three partons at LO. The response of the detector to particles was modeled either with a full ATLAS detector simulation [72] using G 4 [73] or with a fast simulation [74]. The fast simulation is based on a parameterization of the performance of the electromagnetic and hadronic calorimeters and on G 4 elsewhere. All background (signal) samples were prepared using the full (fast) detector simulation. All simulated events were generated with a varying number of minimum-bias interactions overlaid on the hard-scattering event to model the multiple proton-proton interactions in the same and nearby bunch crossings. The minimumbias interactions were simulated with the soft QCD processes of P 8.186 using the A2 tune [75] and the MSTW2008LO PDF set [76]. Corrections were applied to the samples to account for differences between data and simulation for trigger, identification and reconstruction efficiencies.
The proton-proton data analyzed in this paper were collected by ATLAS during 2015 and 2016 at a center-of-mass energy of 13 TeV with up to 50 simultaneous interactions per proton bunch crossing. After application of data-quality requirements related to the beam and detector conditions, the total integrated luminosity corresponds to 36.1 fb −1 . The uncertainty in the combined 2015 and 2016 integrated luminosity is 3.2%. It is derived from a calibration of the luminosity scale using x-y beam-separation scans. This methodology is further detailed in Ref. [77].
The data were collected using the higher-level triggers that select events based on the magnitude of the missing transverse momentum, E miss T . The triggers used are close to fully efficient for events with an offline-reconstructed E miss T greater than 200 GeV.

Event reconstruction
In each event, proton-proton interaction vertices are reconstructed from at least two tracks, each with a transverse momentum p T > 400 MeV and consistent with the beamspot envelope. The primary vertex (PV) of the event is selected as the vertex with the largest p 2 T of the associated tracks. A distinction is made between preselected and signal leptons and jets. Preselected leptons and jets are used in the E miss T computation and are subject to a series of basic quality requirements. Signal leptons and jets are a subset of the preselected objects with more stringent requirements and are used for the definition of signal, control and validation regions. To avoid double-counting of the preselected jets, electrons, and muons, a sequence of overlap-removal procedures based on the angular distance ∆R = (∆y) 2 + (∆φ) 2 is applied. First, any jet reconstructed within ∆R < 0.2 of a preselected electron is rejected. This prevents electromagnetic energy clusters simultaneously reconstructed as an electron and a jet from being selected twice. Next, to remove bremsstrahlung from muons followed by a photon conversion into electron pairs, electrons within ∆R < 0.01 from a preselected muon are discarded. Subsequently, the contamination from muons from decays of heavy hadrons is suppressed by removing muons that are within ∆R < min(0.04 + (10 GeV)/p µ T , 0.4) from preselected jets meeting the previous criteria, or ∆R < 0.2 from a b-tagged jet or a jet containing more than three tracks with p T > 500 MeV. In the former case, the p T -dependent angular separation mitigates the rejection of energetic muons close to jets in boosted event topologies. Finally, jets reconstructed with ∆R < 0.2 from a preselected muon are rejected.
Signal electrons are required to satisfy the likelihood-based tight identification criteria detailed in Ref. [89]. Signal muons and electrons satisfy a sequence of ηand p T -dependent isolation requirements on trackingbased and calorimeter-based variables, defined as the GradientLoose [90] isolation criteria. Compatibility of the signal lepton tracks with the PV is enforced by requiring the distance |z 0 sin θ| to be less than 0.5 mm, where z 0 is the longitudinal impact parameter. In addition, the transverse impact parameter, d 0 , divided by its uncertainty, σ(d 0 ), must satisfy |d 0 /σ(d 0 )| < 3 for signal muons and |d 0 /σ(d 0 )| < 5 for signal electrons.
Corrections derived from data control samples are applied to simulated events to calibrate the reconstruction and identification efficiencies, the momentum scale and resolution of leptons and the efficiency and mis-tag rate of b-tagged jets.

Event selection
Each event must satisfy the trigger selection criteria, and must contain a reconstructed primary vertex. Non-collision background and detector noise are suppressed by rejecting events with any preselected jet not satisfying a set of quality criteria [91]. Exactly one signal lepton, either an electron or a muon, is required. Events with additional preselected leptons are rejected to suppress the dilepton tt, single-top (Wt-channel), Z+jets and diboson backgrounds. The following observables are used in the definition of signal regions in the analysis.
The missing transverse momentum, E miss T , is defined as the magnitude of p miss T , the negative vectorial sum of the transverse momenta of preselected muons, electrons, jets, and identified and calibrated photons. The calculation of p miss T also includes the transverse momenta of all tracks originating from the PV and not associated with any identified object [92, 93].
The transverse mass, m T , is defined from the lepton transverse momentum p T and p miss where ∆φ( p T , p miss T ) is the azimuthal angle between p T and p miss T . For W+jets and semileptonic tt events, in which one on-shell W boson decays leptonically, the observable has an upper endpoint at the W-boson mass. The m T distribution for signal events extends significantly beyond the distributions of the W+jets and semileptonic tt events.
The effective mass, m eff , is the scalar sum of the p T of the signal lepton and all signal jets and E miss T : The effective mass provides good discrimination against SM backgrounds, especially for the signal scenarios where energetic jets are expected. Gluino production leads to higher jet multiplicity than squark production. High-mass sparticles tend to produce harder jets than low-mass sparticles. Thus the optimal m eff value depends on the different signal scenarios. To achieve sensitivity to a wide range of SUSY scenarios with a limited number of signal regions, this variable is binned in the final region definition instead of one simple m eff cut. The detailed description can be found in Section 5.1.
The transverse momentum scalar sum, H T , is defined as where the index j runs over all the signal jets in the event. Empirically, the experimental resolution of E miss T scales with H T , and the ratio E miss T / H T is useful for suppressing background events with large E miss T due to jet mismeasurement.
The aplanarity is a variable designed to provide more global information about the full momentum tensor of the event. It is defined as (3/2) × λ 3 , where λ 3 is the smallest eigenvalue of the normalized momentum tensor [94] calculated using the momenta of the jets and leptons in the event. Typical measured aplanarity values lie in the range 0 -0.3, with values near zero indicating relatively planar background-like events. Signal events tend to have high aplanarity values, since they are more spherical than background events due to multiple objects emitted in the sparticles decay chains.

Signal region definitions
Five sets of event selection criteria, each defining a signal region (SR), are designed to maximize the signal sensitivity. Each SR is labeled by the minimum required number of jets and, optionally, the characteristics of the targeted supersymmetric mass spectrum. Four of the five SRs, 2J, 4J high-x, 4J low-x, and 6J, target the gluino/squark one-step models. The fifth SR, 9J, targets the gluino two-step and pMSSM models. Table 2 summarizes the four SRs targeting the gluino/squark one-step models. The four SRs are mutually exclusive. For setting model-dependent exclusion limits ("excl"), each of the four SRs is further binned in b-veto/b-tag and m eff , and a simultaneous fit is performed across all 28 bins of the four SRs. This choice enhances the sensitivity to a range of new-physics scenarios with different properties such as the presence or absence in the final state of jets containing b-hadrons, and different mass separations between the supersymmetric particles. For model-independent limits and null-hypothesis tests ("disc" for discovery), the event yield above a minimum value of m eff in each SR is used to search for an excess over the SM background.
The 2J SR provides sensitivity to scenarios characterized by a relatively heavyχ 0 1 and small differences between mg, mχ ± 1 , and mχ 0 1 , where most of the decay products tend to have small p T . Events with one low-p T lepton and at least two jets are selected. The minimum lepton p T is 7 (6) GeV for the electron (muon), and the maximum p T is scaled with the number of signal jets in the event as 5 GeV × N jet up to 35 GeV. The maximum p T requirement balances background rejection and signal acceptance for models with increasing mass splittings, where there are more energetic lepton and jets. Stringent requirements on E miss T and on m eff enhance the signal sensitivity by selecting signal events in which the final-state neutralinos are boosted against energetic initial-state radiation (ISR) jets. The SM background is further suppressed by a tight requirement on E miss T /m eff . The 4J high-x SR is optimized for models where mχ 0 1 is fixed to 60 GeV and x ≈ 1, i.e., mχ ± 1 is close to mg. The W boson produced in the chargino decay is significantly boosted, giving rise to a high-p T lepton. The main characteristics of signal events in this model are large m T values and relatively soft jets emitted from the sparticle decay. Tight requirements are placed on E miss T , m T , and E miss T /m eff . The 4J low-x SR targets models where mχ 0 1 is fixed to 60 GeV and x ≈ 0, i.e., mχ ± 1 is close to mχ 0 1 . The large mg /q -mχ ± 1 mass splitting leads to high jet activity, where events are expected to have higher m eff and larger aplanarity than in the high-x scenarios. The W boson tends to be off-shell, leading to small m T , and accordingly an upper bound is imposed to keep this region orthogonal to the 4J high-x SR.  for gluino (squark) Table 3: Overview of the selection criteria for the signal region used for pMSSM and gluino two-step models. SR 9J The 6J SR is optimized for models with x = 1/2, targeting scenarios with large sparticle mass. Events with one high-p T lepton and six or more jets are selected. Requirements on m T , E miss T , m eff , and aplanarity are imposed to reduce the SM background from tt and W + jets production. The sensitivity is improved for scenarios with large mg /q and small mχ 0 1 by introducing a higher m eff bin. Finally, one signal region, 9J SR, is defined to target the pMSSM and gluino two-step models. The selection criteria are summarized in Table 3. At least nine jets are required, targeting the models' long decay chains in which multiple vector or Higgs bosons are produced. The background is further suppressed by tight requirements on the aplanarity and on E miss T / H T . For setting model-dependent exclusion limits ("excl"), the SR is separated into 1000 < m eff < 1500 GeV and m eff > 1500 GeV to achieve good discrimination power for different gluino masses. For model-independent null-hypothesis tests ("disc"), events selected with m eff > 1500 GeV are used to search for an excess over the SM background.

Background estimation
The dominant SM backgrounds in most signal regions originate from top quark (tt and single top) and W+jets production. In this section, the techniques employed to estimate the contribution of these backgrounds in the signal regions are detailed.
Additional sources of background in all signal regions originate from the production of Z+jets, tt in association with a W or Z boson, and diboson (WW, W Z, Z Z) events. Their contributions are estimated entirely using simulated event samples normalized to NLO cross-sections.
The contribution from multi-jet processes with a misidentified lepton is found to be negligible once the lepton isolation and E miss T requirements used in this search are imposed. A data-driven matrix method, following the implementation described in Ref. [21], determined this in previous iterations of the analysis [22]. As this background is found to be negligible, it is not further considered in the analysis.
The dominant top quark and W+jets backgrounds in the 2J, 4J high-x, 4J low-x, and 6J signal regions are estimated by simultaneously normalizing the predicted event yields from simulation to the number of data events observed in dedicated control regions (CR) using the fitting procedure described in Section 8. The simulation is then used to extrapolate the measured background rates to the corresponding signal regions.
The CRs are designed to have high purity in the background process of interest, a sufficiently large number of events to obtain small statistical uncertainties in the background prediction, and a small contamination by events from the signal models under consideration. Moreover, they are designed to have kinematic properties resembling as closely as possible those of the signal regions, in order to provide good estimates of the kinematics of background processes there. This procedure limits the impact of potentially large systematic uncertainties in the expected yields from the extrapolation.
Tables 4-7 list the criteria that define the control regions corresponding to signal regions 2J, 4J highx, 4J low-x, and 6J. As described in Section 5, these signal regions contain multiple bins in m eff . The same binning is maintained for the control regions, so that every signal region bin in m eff has corresponding control regions with the same requirements on m eff and, therefore, the backgrounds are estimated independently in each m eff bin.
Dedicated top and W+jets control regions, respectively denoted by TR and WR, are constructed in each bin of m eff . The TR and WR are distinguished by requiring at least one or exactly zero b-tagged signal jets, respectively. Cross-contamination from top and W+jets processes between these two types of control regions is accounted for in the fit. The measured top and W+jets background rates from the TR and WR regions in a given m eff bin are extrapolated to the signal region within the same m eff bin. The signal regions in a given m eff bin may be further separated into regions with at least one or exactly zero b-tagged signal jets as described in Section 5. For such signal regions separated by b-tagged jet multiplicity, the extrapolation is performed from both the TR and WR regions to each individual bin of b-tagged jet multiplicity.
To validate the extrapolation from control to signal regions using simulated event samples, dedicated validation regions (VRs) are defined for each set of control and signal regions. The selection criteria defining these VRs are also shown in Tables 4-7. The same binning in m eff used in the control and signal regions is also maintained in the validation regions. The VRs are designed to be kinematically close to the signal regions, with only a small contamination from the signal in the models considered in this search. The VRs are not used to constrain parameters in the fit, but provide a statistically independent cross-check Table 4: Overview of the control and validation region selection criteria corresponding to the 2J SR. The top and W+jets control regions are denoted by TR and WR, respectively.

2J
WR   of the extrapolation. The observed event yields in the VRs are found to be consistent with the background prediction as further discussed in Section 8.
One of the dominant background components in the 2J, 4J high-x, 4J low-x, and 6J SRs is tt production with dileptonic final state, where one lepton fails to be reconstructed ("missing lepton") or is a semihadronically decaying τ lepton; this background is characterized by high values of m T . To validate the above described background estimation technique, which is largely a simulation-based extrapolation from low-m T control regions populated by events with semileptonic tt decays, an alternative method was developed. This method (hereafter referred to as the object replacement method) uses events in a dileptonic control region. To emulate the missing lepton case, the p T of one of the two leptons is added vectorially to the calculation of E miss T . To emulate the hadronic τ decay case, one of the two leptons is re-simulated as a hadronic tau decay using the Tauola generator [95] with appropriate energy scale and resolution corrections. The accuracy of this alternative background estimation technique was validated on simulated samples as well as in data validation regions. The background estimates derived from this object replacement method are found to be consistent with those obtained from the standard semi-data-driven approach as further demonstrated in Section 8.
While the background estimation strategy described above works well for the signal regions 2J, 4J    high-x, 4J low-x, and 6J, it is not viable for the 9J SR. The reason for this is that the simulation-based extrapolation from the control regions, which are typically located around the peak region of the transverse mass distribution (m T ∼ 80 GeV), to the high-m T signal regions (m T 80 GeV) is affected by large theoretical uncertainties at high jet multiplicities. Because the peak and tail regions of the m T distribution are dominated by semileptonic and dileptonic final states from tt decays, respectively, additional jets from initial-or final-state radiation are required to obtain the same jet multiplicity for dileptonic tt final states. Inadequate modeling of such additional jets is the dominant source of the theoretical uncertainty. To reduce the dependence on the modeling of additional jets, a dedicated data-driven background estimation technique was designed for the 9J SR. The method relies on the assumption that the m T distribution is approximately invariant under changes in the jet multiplicity requirements. This assumption is found to be valid when tight m eff requirements as used in this analysis are applied such that the overall activity in the calorimeter and thus the missing transverse momentum resolution are not significantly affected by variations in the jet multiplicity. Based on the m T invariance, mutually exclusive control regions CR A,B,C are defined in the m T -N jet plane, where CR A is located at high m T and low N jet , CR B at low m T and low N jet , and CR C at low m T and high N jet . The precise requirements of these regions are defined in Table 8 and illustrated in Figure 2. Based on these regions, the background in the high m T and high N jet signal region can then be estimated with the following equation Table 8: Overview of the control and validation region selection criteria corresponding to the 9J SR. The control regions CR A,A ,B,C,C are further divided into bins of exactly 0 or ≥ 1 b-tagged signal jets to enrich top and W+jets backgrounds, respectively.
9J   Table 3 (signal region) and where N est <region> is the (estimated) number of events in a given region. The residual small correlations between m T and N jet that bias the background estimate in the signal region can then be expressed in terms of a simulation-based closure parameter defined as where N sim <region> is the number of events in a given region as predicted by simulation while N sim,est SR 9J is the estimated number of events in the signal region based on the simulation predictions in regions A, B, and C. The estimated number of background events in the signal region can then be rewritten as where N obs <region> is the observed number of events in a given region, µ C is the normalization parameter in region C, and the normalization parameter µ A/B is fitted simultaneously with the normalization µ B of the backgrounds in region CR B according to The control regions listed in Table 8 are optimized to provide a sufficient number of events in the backgrounds of interest, low contamination from the signal models considered, and a closure parameter f closure close to unity. All control regions are fitted simultaneously in two bins requiring either zero or at least one b-tagged signal jet to enrich the contributions from the W+jets and top backgrounds, respectively. Therefore, the normalization factors µ B , µ C , and µ A/B exist separately for the W+jets and top backgrounds. The top backgrounds considered in the fit comprise tt as well as single-top production processes, which are treated with a common set of normalization parameters.
To validate that the fitted ratio of low-m T to high-m T events (µ A/B ) extrapolates to high values of N jet , a validation region VR m T with seven or eight jets and high m T requirements is introduced. Similarly, a validation region VR N jet with at least nine jets and moderate m T requirements is introduced to validate the extrapolation of the normalization factor µ C in region CR C to higher m T values. Since the normalization factors for different jet multiplicities are expected to differ, a control region CR C along with its normalization factor (µ C ) is introduced. This region is only used to obtain the background estimate in VR m T . Similarly, a control region CR A is constructed to obtain the normalization factor µ A /B that is needed for the background estimation in validation region VR N jet . The definition of the validation regions along with their corresponding control regions is given in Table 8.

Systematic uncertainties
Experimental and theoretical sources of systematic uncertainty are described in this section. Their effects are evaluated for all simulated signal and background events.
The dominant experimental systematic effects are the uncertainties associated with the jet energy scale (JES) and resolution (JER) and with the b-tagging efficiency and mis-tagging rate. The impact of the jet-related uncertainties on the total background prediction ranges from 1.3% in the 6J SR to 18% in the 9J SR. Similarly, the impact of the uncertainties associated with the b-tagging procedure amounts to 1.9% in the 6J SR bins with at least one b-tagged jet and increases to 9.5% in the 6J SR bins with no b-tagged jets. The simulation is reweighted to match the distribution of the average number of proton-proton interactions per bunch crossing (µ) observed in data. The uncertainty in µ is propagated by varying up and down the reweighting factor: it becomes relevant in the signal regions characterized by the highest jet multiplicities.
Uncertainties in the theoretical predictions and the modeling of simulated events are also considered. For the W+jets and the tt and single top backgrounds, they affect the extrapolation from each m eff bin in the control regions to the corresponding bin in the signal regions. In the 9J SR the f closure parameter used in the background estimation in this channel is affected as well. For all the other background sources, they impact the inclusive cross-section of each specific process, the acceptance of the analysis selection requirements and the shape of the m eff distribution in each SR.
An uncertainty stems from the choice of MC event generator modeling the tt, single top, diboson and W/Z+jets processes. For tt and single top, P -B is compared with MG5_ M @NLO [38] and the relative difference in the extrapolation factors is evaluated. For W/Z+jets, the predictions from S are compared with MG5_ M @NLO [38]. For dibosons, the event yield predictions from S are compared with -interfaced to P . The impact of varying the amount of initialand final-state radiation is evaluated for tt and single top production. Specific samples are used, with altered renormalization and factorization scales as well as parton shower and NLO radiation [50]. Moreover, the difference between the predictions from -interfaced to P and to H ++ [96] is computed to estimate the uncertainty associated with the parton shower modeling. For W/Z+jets samples, the uncertainties in the renormalization, factorization, resummation scales and the matching scale between matrix elements and parton shower (CKKW-L) are evaluated by varying up and down by a factor of two the corresponding parameters in S . For tt and W+jets samples, the uncertainties due to choosing the PDF set CT10 [52] are considered.
Inclusive WW bb events generated using MG5_ MC@NLO [38] are compared to the sum of tt and Wt production, to assign an uncertainty to the interference effects between single top and tt production at NLO. The uncertainty in the inclusive Z+jets cross-section, amounting to 5%, is accounted for [97]. An overall 6% systematic uncertainty in the inclusive cross-section of diboson processes is also considered. In Table 9: Breakdown of the dominant systematic uncertainties in the background estimates in the 2J and 4J highx SRs. The individual uncertainties can be correlated and do not necessarily add up in quadrature to the total background uncertainty. The percentages show the size of the uncertainty relative to the total expected background.   addition, the S parameters controlling the renormalization, factorization, resummation and matching scales are varied by a factor of two to estimate the corresponding uncertainties. An uncertainty of 30% is assigned to the small contributions of tt + W/Z/WW.

Signal region
The total systematic uncertainty in the predicted background yields in the various signal regions ranges from 12% in the 2J SR bins with ≥ 1 b-tagged jet, to 50% in the 9J SR. The largest uncertainties in the SR bins with ≥ 1 b-tagged jet originate from the modeling of tt events and amount to 5% in the 2J SR, increasing to 40% in the 9J SR. Similarly, in the SR bins where b-tagged jets are vetoed, the dominant source of systematic uncertainty is the modeling of W+jets events, ranging from 9% in the 6J SR to 20% in the 4J low-x SR. Other important uncertainties are those associated with the finite size of the MC samples, which amount to 18% in the 6J SR, and the theoretical uncertainties originating from the modeling of the diboson background, amounting to 26% in the 6J SR. Tables 9-11 list the breakdown of the dominant systematic uncertainties in background estimates in the various signal regions.
For the signal processes, the modeling of initial-state radiation can be affected by sizable theoretical uncertainty. The uncertainties in the expected yields for SUSY signal models are estimated with variations of a factor of two to the MG5_ MC@NLO parameters corresponding to the renormalization, factorization Table 11: Breakdown of the dominant systematic uncertainties in the background estimates in the 9J SR. The individual uncertainties can be correlated and do not necessarily add up in quadrature to the total background uncertainty. The percentage shows the size of the uncertainty relative to the total expected background.

Signal region 9J
Total background expectation 7 Total background systematic uncertainty ±4 [50%] Theoretical uncertainty ±4 Normalization uncertainty ±2.0 Experimental uncertainty ±1.9 Statistical uncertainty of MC samples ±0.7 and jet matching scales, and to the P shower tune parameters. The overall uncertainties range from about 1% for signal models with large mass splitting between the gluino or squark, the chargino, and the neutralino, to 35% for models with very compressed mass spectra.

Results and Interpretation
The statistical interpretation of the results is performed based on a profile likelihood method [98] using the HistFitter framework [99]. The likelihood function consists of a product of Poisson probability density functions for the signal and control regions that contribute to the fit. The inputs to the likelihood function are the observed numbers of data events and the expected numbers of signal and SM background events in each region. Three normalization factors, one for signal, one for W + jets, one for tt and single top, are introduced to adjust the relative contributions of the main background and signal components. The small sources of SM background, i.e., diboson, Z + jets and tt + V, are estimated directly from simulation. The uncertainties are implemented in the fit as nuisance parameters, which are correlated between the SRs and the CRs. The systematic uncertainties described in Section 7 are constrained by Gaussian probability density functions, while the statistical uncertainties are constrained by Poisson probability density functions.
The observed numbers of events in the signal regions are given in Tables 12-14, along with the SM background prediction as determined with the background-only fit. In a background-only fit, the data event yields in the CRs are used to determine the two background normalization factors: for W + jets and for tt and single top production. The fit is independent of the observation in the SR, and does not consider signal contamination in the CRs. The above-mentioned signal normalization parameter is therefore not included in this fit configuration.
The compatibility of the observed and expected event yields in both the validation and signal regions is illustrated in Figures 3-7. No significant excess in data is observed over the SM prediction.
The top and W + jets background normalization factors obtained for the 2J, 4J low-x, 4J high-x, and 6J SRs are shown in bins of m eff in Figure 8. A trend toward smaller normalization factors at large values of m eff is observed, which demonstrates the necessity of applying the same binning requirements in control and  [GeV]. Uncertainties in the fitted background estimates combine statistical (in the simulated event yields) and systematic uncertainties. The uncertainties in this table are symmetrized for propagation purposes but truncated at zero to remain within the physical boundaries.

2J b-tag
All   signal regions. The predicted event yields from tt events in which both top quarks decay semileptonically are cross-checked using the alternative object-replacement method described in Section 6. Figure 9 shows that the background estimates obtained from the two methods are consistent. Figures 10-11 show the m eff distribution in 2J, 4J low-x, 4J high-x and 6J in b-tag and b-veto signal regions after fit. Figure 12 shows the m eff distribution in 9J signal region after fit. The uncertainty bands plotted include all statistical and systematic uncertainties. The dashed lines stand for the benchmark signal samples.
Using the results of the background-only fit, a model-independent limit fit is performed to test for the presence of any beyond-the-Standard-Model (BSM) physics processes that contribute to the SR ("disc" SR in Table 2) . The BSM signal is assumed to contribute only to the SR and not to the CRs, thus giving a conservative estimate of background in the SR. Observed (S 95 obs ) and expected (S 95 exp ) 95% confidence level (CL) upper limits on the number of BSM signal events are derived using the CL s prescription [100]. Table 15 presents these limits, together with the upper limits on the visible BSM cross-section, σ 95 obs , defined as the product of acceptance, selection efficiency and production cross-section. The upper limits on the visible BSM cross-section are calculated by dividing the observed upper limit on the beyond-SM events by the integrated luminosity of 36.1 fb −1 . Moreover, the discovery p-values are given. They quantify the probability under the background-only hypothesis to produce event yields greater than or equal to the observed data.
Additionally, the results are interpreted in the specific supersymmetric scenarios described in Section 3 using model-dependent limit fits. A model-dependent limit fit takes the data event yields in multiple, statistically independent SRs and their associated CRs to compute an upper limit on the cross-section of a targeted SUSY model. The fit includes the expected signal contributions to the SRs and to the CRs, scaled by a floating signal normalization factor. The background normalization factors are also determined simultaneously in the fit.
The sparticle mass in a specific SUSY model can be excluded if the upper limit of the signal normalization factor obtained in the fit is smaller than unity.
For the gluino/squark one-step models, a model-dependent fit is performed over all bins of the 2J, 4J high-x, 4J low-x, and 6J SRs. An independent set of background normalization factors are allocated for each bin of each SR ("excl" SR in Table 2) and its associated CRs. Figure 13 (top and middle) shows the observed and expected exclusion bounds at 95% CL for the one-step simplified models with gluino and squark production. Gluino masses up to 2.1 TeV and squark masses up to 1.25 TeV are excluded. Figure 13 (bottom) shows the exclusion contours of the 9J SR (Table 3) for the gluino two-step as well as the pMSSM scenario described in Section 3. In both cases the limits reach well beyond 1.7 TeV in gluino mass.  Table 15: Results of the model-independent limit fits. For each SR, the observed 95% CL upper limit on the visible cross-section ( σ 95 obs ), the observed (S 95 obs ) and expected (S 95 exp ) 95% CL upper limits on the BSM event yield, and the one-sided discovery p-value (p(s = 0)) are presented. The p-values are capped at 0.5 if fewer events than the fitted background estimate are observed.

SR
V R a p la n a ri ty V R a p la n a ri ty m V R a p la n a ri ty m V R a p la n a ri ty m V R a p la n a ri ty V R a p la n a ri ty m V R a p la n a ri ty m V R a p la n a ri ty m V R a p la n a ri ty m        Figure 13: Exclusion contours for gluino one-step x = 1/2 (top left), gluino one-step variable-x (top right), squark one-step x = 1/2 (middle left) and squark one-step variable-x (middle right), gluino two-step (bottom left), and the pMSSM scenario (bottom right). The red solid line corresponds to the observed limit with the red dotted lines indicating the ±1σ variation of this limit due to the effect of theoretical scale and PDF uncertainties in the signal cross-section. The dark gray dashed line indicates the expected limit with the yellow band representing the ± 1 σ variation of the median expected limit due to the experimental and theoretical uncertainties. For reference, exclusion bounds from previous searches with 20.3 fb −1 at 8 TeV center-of-mass energy [28] and 3.2 fb −1 at 13 TeV center-ofmass energy [22,30] are overlaid where applicable by the gray area (the observed limit is shown by the solid line, while the dashed line shows the expected limit) .

Conclusion
A search for the pair production of squarks and gluinos in proton-proton collisions provided by the LHC at a center-of-mass energy of √ s = 13 TeV has been performed by the ATLAS Collaboration. Events containing one isolated electron or muon, two or more jets, and large missing transverse momentum are selected in the data collected in 2015 and 2016, corresponding to an integrated luminosity of 36.    [79] M. Cacciari, G. P. Salam      [94] C. Chen, New approach to identifying boosted hadronically decaying particles using jet substracture in its center-of-mass frame, Phys.