1 Introduction

In the standard model (SM), the production of multiple jets in conjunction with two same-sign (SS) or three or more charged leptons is a very rare process in proton-proton (\({\text {p}}{\text {p}}\)) collisions. These final states provide a promising starting point in the search for physics beyond the SM (BSM). Many models attempting to address the shortcomings of the SM lead to such signatures. Examples include the production of supersymmetric (SUSY) particles [1, 2], SS top quark pairs [3, 4], scalar gluons (sgluons) [5, 6], heavy scalar bosons of extended Higgs sectors [7, 8], Majorana neutrinos [9], and vector-like quarks [10].

In SUSY models [11,12,13,14,15,16,17,18,19], the decay chain of pair-produced gluinos or squarks can contain multiple \({\text {W}}\) or \({\text {Z}}\) bosons, with the potential to have at least one pair of SS \({\text {W}}\) bosons. Such a decay chain is realized, for example in gluino pair production, when a gluino decays into a top quark-antiquark pair and a neutralino, or into a pair of quarks and a chargino that subsequently decays into a \({\text {W}}\) boson and a neutralino. In R parity [20] conserving (RPC) scenarios, the lightest SUSY particle is neutral and stable and escapes detection, leading to an imbalance in the measured transverse momentum. The magnitude of the missing transverse momentum strongly depends on the details of the model, and in particular on the mass spectrum of the particles involved. Scenarios with R parity violation (RPV) [21, 22] additionally allow decays of SUSY particles into SM particles only, leading in many cases to signatures with little or no missing transverse momentum. For many SUSY models, the SS and multilepton signatures provide complementarity with searches in the zero- or one-lepton final states, and they are particularly suitable for probing compressed mass spectra and other scenarios involving low-momentum leptons or low missing transverse momentum. Both the ATLAS [23] and CMS [24, 25] Collaborations have carried out searches in these channels using LHC data collected up to and including 2016. The ATLAS Collaboration has also recently released a search with the full data set recorded between 2015 and 2018 [26].

In this paper, we extend and refine the searches described in Refs. [24, 25] using a larger data set of \({\text {p}}{\text {p}}\) collisions at \(\sqrt{s} = 13\,{\text {TeV}} \) recorded by the CMS detector at the CERN LHC in 2016–2018, corresponding to an integrated luminosity of \(137{\,{\text {fb}}^{-1}} \). We base our search on an initial selection of events with at least two hadronic jets and two SS or three or more light leptons (electrons and muons), including those from leptonic decays of \(\tau \) leptons. Several signal regions (SRs) are then constructed with requirements on variables such as the number of leptons, the number of jets (possibly identified as originating from \({\text {b}}\) quarks), and the magnitude of missing transverse momentum. A simultaneous comparison of the observed and SM plus BSM expected event yields in all SRs is performed to constrain the BSM models described in Sect. 2. After a brief description of the CMS experiment in Sect. 3, we present the details of the search strategy and event selection in Sect. 4 and discuss the various relevant backgrounds from SM processes in Sect. 5. The systematic uncertainties considered in the analysis are presented in Sect. 6. In Sect. 7, the observed yields are compared to the background expectation and the results are interpreted to constrain the various BSM models introduced earlier. Model independent limits are also derived. Finally, the main results are summarized in Sect. 8.

2 Background and signal simulation

Monte Carlo (MC) simulations are used to study the SM backgrounds and to estimate the event selection efficiency of the BSM signals under consideration. Three sets of simulated events for each process are used in order to match the different data taking conditions in 2016, 2017, and 2018.

The hard scattering process of the dominant backgrounds estimated from simulation (including the \(\hbox {t}{\bar{\hbox {t}}} {\text {W}}\), \(\hbox {t}{\bar{\hbox {t}}} {\text {Z}}\) and \({\text {W}}{\text {Z}}\) contributions) is simulated with the MadGraph 5_amc@nlo 2.2.2 (2.4.2) [27,28,29] generator for 2016 (2017 and 2018) conditions. An exception is the \({\text {W}}{\text {Z}}\) process for the 2016 conditions that, as with a few subdominant backgrounds, is simulated using the powheg  v2 [30,31,32,33,34] next-to-leading order (NLO) generator. Samples of signal events, as well as of SS \({\text {W}}\) boson pairs and other very rare SM processes, are generated at leading order (LO) accuracy with MadGraph 5_amc@nlo, with up to two additional partons in the matrix element calculations. The set of parton distribution functions (PDFs) used was NNPDF3.0 [35] for the 2016 simulation and NNPDF3.1 [36] for the 2017 and 2018 simulations.

Parton showering and hadronization, as well as the double parton scattering production of \({\text {W}}^{\pm } {\text {W}}^{\pm }\), are described using the pythia  8.230 generator [37] with the CUETP8M1 (CP5) underlying event tune for 2016 (2017 and 2018) simulation [38,39,40]. The response of the CMS detector is modeled using the Geant4 program [41] for SM background samples, while the CMS fast simulation package [42, 43] is used for signal samples.

To improve the MadGraph 5_amc@nlo modeling of the multiplicity of additional jets from initial-state radiation (ISR), 2016 MC events are reweighted according to the number of ISR jets (\(N_J^\text {ISR}\)). The reweighting factors are extracted from a study of the light-flavor jet multiplicity in dilepton \(\hbox {t}{\bar{\hbox {t}}}\) events. They vary between 0.92 and 0.77 for \(N_J^\text {ISR}\) between 1 and 4, with one half of the deviation from unity taken as the systematic uncertainty. This reweighting is not necessary for the 2017 and 2018 MC samples that are generated using an updated pythia tune.

The phenomenology of a given SUSY model strongly depends on its underlying details such as the masses of the SUSY particles and their couplings with the SM particles and each other, many of which can be free parameters. The signal models used by this search are simplified SUSY models [44, 45] of either gluino or squark pair production, followed by a variety of RPC (Figs. 1, 2) or RPV (Fig. 3) decays and where several leptons can arise in the final state. Production cross sections are calculated at approximate next-to-next-to-leading order plus next-to-next-to-leading logarithmic (NNLO+NNLL) accuracy [46,47,48,49,50,51,52,53,54,55,56,57,58]. The branching fractions for the decays shown are assumed to be 100%, unless otherwise specified, and all decays are assumed to be prompt.

Gluino pair production models giving rise to signatures with up to four \({\text {b}}\) quarks and up to four \({\text {W}}\) bosons are shown in Fig. 1. In these models, the gluino decays to the lightest squark (\({\tilde{\text {g}}} \rightarrow {\tilde{\text {q}}} {\text {q}}\)), which in turn decays to same-flavor (\({\tilde{\text {q}}} \rightarrow {\text {q}}{\tilde{\chi }}^{0}_{1} \)) or different-flavor (\({\tilde{\text {q}}} \rightarrow {\text {q}}' {\tilde{\chi }}^\pm _{1} \)) quarks. The chargino (\({\tilde{\chi }}^\pm _{1} \)) decays to a \({\text {W}}\) boson and a neutralino (\({\tilde{\chi }}^{0}_{1} \)) via \({\tilde{\chi }}^\pm _{1} \rightarrow {\text {W}}^{\pm } {\tilde{\chi }}^{0}_{1} \), where the \({\tilde{\chi }}^{0}_{1}\) is taken to be the lightest stable SUSY particle and escapes detection.

The first scenario, denoted by \(\mathrm {T}1{{\text {t t t t}}}\) and displayed in Fig. 1a, includes an off-shell top squark (\({\tilde{\text {t}}} \)) leading to the three-body decay of the gluino, \({\tilde{\text {g}}} \rightarrow \hbox {t}{\bar{\hbox {t}}} {\tilde{\chi }}^{0}_{1} \), resulting in events with four \({\text {W}}\) bosons and four \({\text {b}}\) quarks. Figure 1b presents a similar model (\(\mathrm {T}5{{\text {t t b b W W}}}\)) where the gluino decay results in a chargino that further decays into a neutralino and a \({\text {W}}\) boson. The model shown in Fig. 1c (\(\mathrm {T}5{{\text {t t t t}}}\)) is the same as \(\mathrm {T}1{{\text {t t t t}}}\) except that the intermediate top squark is on-shell. The mass splitting between the \({\tilde{\text {t}}} \) and the \({\tilde{\chi }}^{0}_{1}\) is taken to be \(m_{{\tilde{\text {t}}}} - m_{{\tilde{\chi }}^{0}_{1}} = m_{{\text {t}}}\), where \(m_{{\text {t}}}\) is the top quark mass. This choice maximizes the kinematic differences between this model and \(\mathrm {T}1{{\text {t t t t}}}\), and also corresponds to one of the most challenging regions of parameter space for the observation of the \({\tilde{\text {t}}} \rightarrow {\text {t}}{\tilde{\chi }}^{0}_{1} \) decay since the neutralino is produced at rest in the top squark rest frame. The decay chain of Fig. 1d (\(\mathrm {T}5{{\text {t t c c}}}\)) is identical to that of \(\mathrm {T}5{{\text {t t t t}}}\) except that the \({\tilde{\text {t}}} \) decay involves a \({\text {c}}\)quark. In Fig. 1e, the decay process includes a virtual light-flavor squark, leading to three-body decays of \({\tilde{\text {g}}} \rightarrow {\text {q q}}' {\tilde{\chi }}^\pm _{1} \) or \({\tilde{\text {g}}} \rightarrow {\text {q q}}' {\tilde{\chi }} _{2}^{0}\), with a resulting signature of two \({\text {W}}\) bosons, two \({\text {Z}}\) bosons, or one of each (the case shown in Fig. 2e), and four light-flavor jets. This model, \(\mathrm {T}5{{\text {qqqqWZ}}}\), with a resulting signature of one \({\text {W}}\) boson and one \({\text {Z}}\) boson, is studied with two different assumptions for the chargino mass: \(m_{{\tilde{\chi }}^\pm _{1}} = 0.5(m_{{\tilde{\text {g}}}} + m_{{\tilde{\chi }}^{0}_{1}})\), and \(m_{{\tilde{\chi }}^\pm _{1}} = m_{{\tilde{\chi }}^{0}_{1}}+20 \,{\text {GeV}} \), producing on- and off-shell bosons, respectively. The model is also considered with the assumption of decays to two \({\text {W}}\) bosons exclusively (\(\mathrm {T}5{{\text {q q q qWW}}}\)).

Figure 2a shows a model of bottom squark production with subsequent decay of \({\tilde{\text {b}}} _{1}\rightarrow {\text {t}}{\tilde{\chi }}^\pm _{1} \), yielding two \({\text {b}}\) quarks and four \({\text {W}}\) bosons. This model, \(\mathrm {T}6{{\text {ttWW}}}\), is considered as a function of the the lightest bottom squark, \({\tilde{\text {b}}} _{1}\), and \({\tilde{\chi }}^\pm _{1}\) masses. The \({\tilde{\chi }}^{0}_{1}\) mass is fixed to be \(50\,{\text {GeV}} \), causing two of the \({\text {W}}\) bosons to be produced off-shell when the \({\tilde{\chi }}^\pm _{1}\) mass is less than approximately \(130\,{\text {GeV}} \). Figure 2b displays a model similar to \(\mathrm {T}6{{\text{ttWW}}}\), but with top squark pair production and a subsequent decay of \({\tilde{\text {t}}} _{2}\rightarrow {\tilde{\text {t}}} _{1}{\text {H}}/{\text {Z}}\), with \({\tilde{\text {t}}} _{1}\rightarrow {\text {t}}{\tilde{\chi }}^{0}_{1} \), producing signatures with two \({\text {H}}\) bosons, two \({\text {Z}}\) bosons, or one of each. In this model, \(\mathrm {T}6{{\text{ttHZ}}}\), the \({\tilde{\chi }}^{0}_{1}\) mass is fixed such that \(m({\tilde{\text {t}}} _{1})-m({\tilde{\chi }}^{0}_{1})=m_{{\text {t}}}\).

The R parity violating decays considered in this analysis are \(\mathrm {T}1{{\text{qqqqL}}} \) (Fig. 3a) and \(\mathrm {T}1{{\text {tbs}}}\) (Fig. 3b). In \(\mathrm {T}1{{\text {qqqqL}}} \), the gluino decays to the lightest squark (\({\tilde{\text {g}}} \rightarrow {\tilde{\text {q}}} {\text {q}}\)), which in turn decays to a quark (\({\tilde{\text {q}}} \rightarrow {\text {q}}{\tilde{\chi }}^{0}_{1} \)), but decays with the \({\tilde{\chi }}^{0}_{1} \) off shell (violating R parity) into two quarks and a charged lepton, giving rise to a prompt 5-body decay of the gluino. In \(\mathrm {T}1{{\text{tbs}}}\), each gluino decays into three different SM quarks (a top, a bottom, and a strange quark).

Fig. 1
figure 1

Diagrams illustrating the simplified RPC SUSY models with gluino production considered in this analysis

Fig. 2
figure 2

Diagrams illustrating the simplified RPC SUSY models with squark production considered in this analysis

Fig. 3
figure 3

Diagrams illustrating the two simplified RPV SUSY models considered in this analysis

3 The CMS detector and event reconstruction

The central feature of the CMS detector is a superconducting solenoid of \(6\,{\text {m}} \) internal diameter, providing a magnetic field of \(3.8\,{\text {T}} \). Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (\(\eta \)) 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. [59].

Events of interest are selected using a two-tiered trigger system [60]. 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\,{\text {kHz}} \) within a fixed time interval of less than 4 μs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around \(1\,{\text {kHz}} \) before data storage.

The reconstructed vertex with the largest value of summed physics-object squared-transverse-momentum is taken to be the primary \({\text {pp}}\) interaction vertex. The physics objects are the jets, clustered using the jet finding algorithm of Refs. [61, 62] with the tracks assigned to the vertex as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the transverse momentum (\(p_{{\mathrm {T}}}\)) of those jets.

The particle-flow (PF) algorithm [63] aims to reconstruct and identify each individual particle in an event, with an optimized combination of information from the various elements of the CMS detector. The energy of photons is directly obtained from the ECAL measurement. The energy of electrons is determined from a combination of the electron momentum at the primary interaction vertex as determined by the tracker, the energy of the corresponding ECAL cluster, and the energy sum of all bremsstrahlung photons spatially compatible with the electron track [64]. The momentum of muons is obtained from the curvature of the corresponding track, combining information from the silicon tracker and the muon system [65]. The energy of charged hadrons is determined from a combination of their momentum measured in the tracker and the matching ECAL and HCAL energy deposits, corrected for the response function of the calorimeters to hadronic showers. The energy of neutral hadrons is obtained from the corresponding corrected ECAL and HCAL energies.

Hadronic jets are clustered from charged PF candidates associated with the primary vertex and from all neutral PF candidates using the anti-\(k_{{\mathrm {T}}}\) algorithm [61, 62] with a distance parameter of 0.4. The jet momentum is determined as the vectorial sum of all PF candidate momenta in the jet. An offset correction is applied to jet energies to take into account the contribution from pileup [66]. Additional jet energy corrections are derived from simulation to bring the detector response to unity, and are improved with in situ measurements of the energy balance in dijet, multijet, photon+jets, and leptonically decaying \({\text {Z+jets}}\) events [67, 68]. Additional selection criteria are applied to each jet to remove jets potentially dominated by instrumental effects or reconstruction failures. Jets originating from \({\text {b}}\) quarks are identified as \({\text {b}}\)-tagged jets using a deep neural network algorithm, DeepCSV [69], with a working point chosen such that the efficiency to identify a \({\text {b}}\) jet is 55–70% for a jet \(p_{{\mathrm {T}}}\) between 20 and \(400\,{\text {GeV}} \). The misidentification rate for a light-flavor jet is 1–2% in the same jet \(p_{{\mathrm {T}}}\) range.

The vector \({\vec {p}}_{{\mathrm {T}}}^{{\text {miss}}}\) is defined as the projection onto the plane perpendicular to the beams of the negative vector sum of the momenta of all reconstructed PF candidates in an event [70]. Its magnitude, called missing transverse momentum, is referred to as \(p_{{\mathrm {T}}} ^{\text {miss}}\). The scalar \(p_{{\mathrm {T}}}\) sum of all jets in an event is referred to as \(H_{\text {T}}\).

4 Search strategy and event selection

The search strategy is similar to the one adopted in Refs. [24, 25]. The event selection requires the presence of at least two hadronic jets and at least two leptons, among which is an SS pair, as described below. Each selected event is assigned to an SR, based on its content. Maximum likelihood fits of the background (or signal plus background) predictions to the data in all SRs are then performed. Such a strategy ensures sensitivity to a broad range of possible signatures of new physics, even beyond the signal benchmarks considered in this analysis.

The kinematic requirements applied to leptons and jets are presented in Table 1. The analysis requires at least two jets with \(p_{{\mathrm {T}}} >40\,{\text {GeV}} \) and two light SS leptons with \(p_{{\mathrm {T}}} >15\,{\text {GeV}} \) (\(10\,{\text {GeV}} \)) for electrons (muons). Electrons are identified based on a discriminant using shower shape and track quality variables, while the muon identification relies on the quality of the geometrical matching between the tracker and muon system measurements. In order to reject leptons from the decay of heavy flavor hadrons, the tracks are required to have an impact parameter compatible with the position of the primary vertex. Several isolation criteria are also applied, based on the scalar sum of hadron and photon \(p_{{\mathrm {T}}}\) within a cone centered on the lepton direction and whose radius decreases with its \(p_{{\mathrm {T}}}\), the ratio of the \(p_{{\mathrm {T}}}\) of the lepton to that of the closest jet, and the relative \(p_{{\mathrm {T}}}\) of the lepton to that of the closest jet after lepton momentum subtraction. These criteria are designed to mitigate the loss of lepton efficiency caused by lepton-jet overlaps that occurs frequently in events with significant hadronic activity. A more detailed description of the set of identification and isolation variables used in the lepton selection can be found in Ref. [71].

The lepton reconstruction and identification efficiency is in the range of 45–70% (70–90%) for electrons (muons), with \(p_{{\mathrm {T}}} >25\,{\text {GeV}} \), increasing as a function of \(p_{{\mathrm {T}}}\) and reaching the maximum value for \(p_{{\mathrm {T}}} >60\,{\text {GeV}} \). In the low-momentum regime, \(15<p_{{\mathrm {T}}} <25\,{\text {GeV}} \) for electrons and \(10<p_{{\mathrm {T}}} <25\,{\text {GeV}} \) for muons, the efficiencies are approximately 40% for electrons and 55% for muons. The lepton trigger efficiency for electrons is in the range of 90–98%, converging to the maximum value for \(p_{{\mathrm {T}}} >30 \,{\text {GeV}} \), and it is around 92% for muons.

Table 1 Transverse momentum and pseudorapidity requirements for leptons and jets. Note that the \(p_{{\mathrm {T}}} \) thresholds to count jets and \({\text {b}}\)-tagged jets are different; the jet multiplicity \(N_\text {jets}\) includes \({\text {b}}\)-tagged jets if their \(p_{{\mathrm {T}}}\) exceeds \(40\,{\text {GeV}} \)

In order to reduce backgrounds from the decays of \({\text {c}}\)- and \({\text {b}}\)-hadrons or from the Drell–Yan process, we reject events with same-flavor lepton pairs with invariant mass (\(m_{\ell \ell }\)) less than \(12\,{\text {GeV}} \), where leptons are reconstructed with a looser set of requirements compared to the nominal selection. Furthermore, events containing a lepton pair with \(m_{\ell \ell }< 8 \,{\text {GeV}} \), regardless of charge or flavor, are rejected in order to emulate a similar condition applied at the trigger level. Events are then separated according to the \(p_{{\mathrm {T}}} \) of the leptons forming the SS pair: high-high if both have \(p_{{\mathrm {T}}} > 25\,{\text {GeV}} \), low-low if both have \(p_{{\mathrm {T}}} < 25\,{\text {GeV}} \), and high-low otherwise.

Two sets of trigger algorithms are used to select the events: pure dilepton triggers, which require the presence of two isolated leptons with \(p_{{\mathrm {T}}} \) thresholds on the leading (subleading) lepton in the 17–23 (8–12) \(\,{\text {GeV}}\) range, and dilepton triggers with no isolation requirements, a lower \(p_{{\mathrm {T}}} \) threshold of \(8\,{\text {GeV}} \), an invariant mass condition \(m_{\ell \ell }> 8 \,{\text {GeV}} \) to reject low mass resonances, and with a minimum \(H_{\text {T}} \) in the range of \(300{-}350\,{\text {GeV}} \). The ranges listed here reflect the varying trigger conditions during the data taking periods. The pure dilepton triggers are used to select high-high and high-low pairs, while low-low pairs are selected using the triggers with \(H_{\text {T}} \) requirements.

Six exclusive categories are then defined as follows:

  • High-High SS pair, significant \(p_{{\mathrm {T}}} ^{\text {miss}}\) (HH): exactly 2 leptons, both with \(p_{{\mathrm {T}}} >25 \,{\text {GeV}} \), and \(p_{{\mathrm {T}}} ^{\text {miss}} >50 \,{\text {GeV}} \);

  • High-Low SS pair, significant \(p_{{\mathrm {T}}} ^{\text {miss}}\) (HL): exactly 2 leptons, one with \(p_{{\mathrm {T}}} >25 \,{\text {GeV}} \), one with \(p_{{\mathrm {T}}} <25 \,{\text {GeV}} \), and \(p_{{\mathrm {T}}} ^{\text {miss}} >50 \,{\text {GeV}} \);

  • Low-Low SS pair, significant \(p_{{\mathrm {T}}} ^{\text {miss}}\) (LL): exactly 2 leptons, both with \(p_{{\mathrm {T}}} <25 \,{\text {GeV}} \) and \(p_{{\mathrm {T}}} ^{\text {miss}} >50 \,{\text {GeV}} \);

  • Low \(p_{{\mathrm {T}}} ^{\text {miss}}\) (LM): exactly 2 leptons, both with \(p_{{\mathrm {T}}} >25 \,{\text {GeV}} \), and \(p_{{\mathrm {T}}} ^{\text {miss}} <50 \,{\text {GeV}} \); and

  • Multilepton with an on-shell Z boson (on-\({\text {Z}}\) ML): \(\ge \)3 leptons, at least one with \(p_{{\mathrm {T}}} >25 \,{\text {GeV}} \), \(p_{{\mathrm {T}}} ^{\text {miss}} >50 \,{\text {GeV}} \), \(\ge {\text {Z}}\) boson candidate formed by a pair of opposite-sign (OS), same-flavor leptons with \(76< m_{\ell \ell }< 106 \,{\text {GeV}} \).

  • Multilepton without an on-shell Z boson (off-\({\text {Z}}\) ML): same as on-\({\text {Z}}\) ML but without a \({\text {Z}}\) boson candidate.

The categories are typically sensitive to different new physics scenarios and enriched in different SM backgrounds. For example the HH category drives the sensitivity for most of the RPC scenarios (\(\mathrm {T}1{{\text {tttt}}}\), \(\mathrm {T}5{{\text {ttbbWW}}}\), \(\mathrm {T}5{{\text {tttt}}}\), \(\mathrm {T}1{{\text {tttt}}}\), \(\mathrm {T}5{{\text {qqqqWW}}}\)) with a large mass splitting. The HL and LL categories become relevant for lower mass splitting when one or both leptons tend to be soft. Scenarios resulting in the presence of one or multiple \({\text {Z}}\) bosons in the final state such as \(\mathrm {T}5{{\text {q}} {\text {qqqWZ}}}\) and \(\mathrm {T}6{{\text {ttHZ}}}\) will typically be primarily constrained by the on-Z or off-Z category, also depending on the considered SUSY mass spectrum. Finally the LM category enhances the analysis sensitivity for RPV scenarios, in particular for \(\mathrm {T}1{{\text {qqqqL}}} \) where no genuine \(p_{{\mathrm {T}}} ^{\text {miss}}\) is expected.

Various SRs are constructed based on the jet multiplicity \(N_\text {jets}\), the \({\text {b}}\)-tagged jet multiplicity \(N_{{\text {b}}}\), \(H_{\text {T}}\), \(p_{{\mathrm {T}}} ^{\text {miss}}\), the charge of the SS pair, and \(m_{\text {T}} ^{\text {min}}\), which is defined below. The \(m_{\text {T}} ^{\text {min}}\) variable, introduced in Ref. [71], is defined as the minimum of the transverse masses calculated from each of the leptons forming the SS pair and \({\vec {p}}_{{\mathrm {T}}}^{{\text {miss}}}\), except for the on-\({\text {Z}}\) ML category where we only consider the transverse mass computed using the leptons not forming the \({\text {Z}}\) candidate. It is characterized by a kinematic cutoff for events where \(p_{{\mathrm {T}}} ^{\text {miss}}\) only arises from the leptonic decay of a single \({\text {W}}\) boson and is effective at discriminating signal and background signatures.

A subset of SRs is split by the charge of the leptons in an SS pair which is used to take advantage of the charge asymmetry in most of the background processes, such as \({\text {W}}{\text {Z}}\), \(\hbox {t}{\bar{\hbox {t}}} {\text {W}}\) or SS \({\text {W}}{\text {W}}\). The SRs corresponding to each category, HH, HL, LL, LM, on-\({\text {Z}}\) ML, and off-\({\text {Z}}\) ML, are summarized in Tables 23456 and 7, respectively. The binning ranges are chosen to maximize the sensitivity to a variety of SUSY benchmark points and are such that the expected SM yield in any SR has relative statistical uncertainties typically smaller than unity.

Table 2 The SR definitions for the HH category. Charge-split regions are indicated with (++) and (- -). The three highest \(H_{\text {T}} \) regions are split only by \(N_\text {jets} \), resulting in 62 regions in total. Quantities are specified in units of \(\,{\text {GeV}}\) where applicable
Table 3 The SR definitions for the HL category. Charge-split regions are indicated with (++) and (- -). There are 43 regions in total. Quantities are specified in units of \(\,{\text {GeV}}\) where applicable
Table 4 The SR definitions for the LL category. All SRs in this category require \(N_\text {jets} \ge 2\). There are 8 regions in total. Quantities are specified in units of \(\,{\text {GeV}}\) where applicable
Table 5 The SR definitions for the LM category. All SRs in this category require \(p_{{\mathrm {T}}} ^{\text {miss}} < 50\,{\text {GeV}} \) and \(H_{\text {T}} >300\,{\text {GeV}} \). The two high-\(H_{\text {T}} \) regions are split only by \(N_\text {jets} \), resulting in 11 regions in total. Quantities are specified in units of \(\,{\text {GeV}}\) where applicable
Table 6 The SR definitions for the on-\({\text {Z}}\) ML category. All SRs in these categories require \(N_\text {jets} \ge 2\). Regions marked with \({}^\dagger \) are split by \(m_{\text {T}} ^{\text {min}} =120\,{\text {GeV}} \), with the high-\(m_{\text {T}} ^{\text {min}} \) region specified by the second SR label. There are 23 regions in total. Quantities are specified in units of \(\,{\text {GeV}}\) where applicable
Table 7 The SR definitions for the off-\({\text {Z}}\) category. All SRs in these categories require \(N_\text {jets} \ge 2\). Regions marked with \({}^\dagger \) are split by \(m_{\text {T}} ^{\text {min}} =120\,{\text {GeV}} \), with the high-\(m_{\text {T}} ^{\text {min}} \) region specified by the second SR label. There are 21 regions in total. Quantities are specified in units of \(\,{\text {GeV}}\) where applicable

5 Backgrounds

Several SM processes can lead to the signatures studied in this analysis. There are three background categories, depending on the lepton content of the event:

  • Events with two or more prompt leptons, including an SS pair;

  • Events with at least one nonprompt lepton (defined below); and

  • Events with a pair of OS leptons, one of which is reconstructed with the wrong charge.

The first category includes a variety of low cross section processes where multiple electroweak bosons are produced, possibly in the decay of top quarks, which then decay leptonically leading to an SS lepton pair. This category usually dominates the background yields in SRs with large \(p_{{\mathrm {T}}} ^{\text {miss}}\) or \(H_{\text {T}}\) and in most of the ML SRs with a \({\text {Z}}\) candidate. The main contributions arise from the production of a \({\text {W}}{\text {Z}}\) or an SS \({\text {W}}\)pair, or of a \(\hbox {t}{\bar{\hbox {t}}} \) pair in association with a \({\text {W}}\), \({\text {Z}}\) or \({\text {H}}\) boson. The event yields for these processes are estimated individually. In contrast, the expected event yields from other rare processes (including \({\text {Z}}{\text {Z}}\), triple boson production, \({\text {tWZ}}\), \({\text {tZq}}\), \(\hbox {t}{\bar{\hbox {t}}} \hbox {t}{\bar{\hbox {t}}} \), and double parton scattering) are summed up into a single contribution denoted as “Rare”. Processes including a genuine photon, such as \({\text {W}}{\upgamma }\), \({\text {Z}}{\upgamma }\), \(\hbox {t}{\bar{\hbox {t}}} {\upgamma }\), and \({\text {t}}{\upgamma }\), are also considered and grouped together. They are referred to as “X\({\upgamma }\)”. All contributions from this category are estimated using simulated samples. Correction factors are applied to take into account small differences between data and simulation, including trigger, lepton selection, and \({\text {b}}\) tagging efficiencies, with associated systematic uncertainties listed in Sect. 6.

The second category consists of events where one of the selected leptons, generically denoted as “nonprompt lepton”, is either a decay product of a heavy flavor hadron or, more rarely, a misidentified hadron. This category is typically the dominant one in SRs with moderate or low \(p_{{\mathrm {T}}} ^{\text {miss}}\) or low \(m_{\text {T}} ^{\text {min}}\) (except for the on-\({\text {Z}}\) ML SRs). This background is estimated directly from data using the “tight-to-loose” method [24, 25]. This method is based on the probability for a nonprompt lepton passing loose selection criteria to also satisfy the tighter lepton selection used in the analysis. The number of events in an SR with N leptons, including at least one nonprompt lepton, can be estimated by applying this probability to a corresponding control region (CR) of events with N loose leptons where at least one of them fails the tight selection.

The measurement of the tight-to-loose ratio is performed in a sample enriched in dijet events with exactly one loose lepton, low \(p_{{\mathrm {T}}} ^{\text {miss}}\), and low \(m_{\text {T}} ^{\text {min}}\). This sample is contaminated by prompt leptons from \({\text {W}}\) boson decays. The contamination is estimated from the \(m_{\text {T}} ^{\text {min}}\) distribution, and it is subtracted before calculating the ratio. The tight-to-loose ratio is computed separately for electrons and muons, and is parameterized as a function of the lepton \(\eta \) and \(p_{{\mathrm {T}}} ^\text {corr}\). The \(p_{{\mathrm {T}}} ^\text {corr}\) variable is defined as the sum of the lepton \(p_{{\mathrm {T}}}\) and the energy in the isolation cone exceeding the isolation threshold value applied to tight leptons. This parametrization improves the stability of the tight-to-loose ratio with respect to variations in the \(p_{{\mathrm {T}}}\) of the partons from which the leptons originate.

The performance of the tight-to-loose ratio was assessed in a MC closure test. A tight-to-loose ratio was extracted from a MC sample of QCD events. This ratio was then used to predict the number of events with one prompt and one nonprompt SS dileptons in MC \(\hbox {t}{\bar{\hbox {t}}}\) and \({\text {W+jets}}\) events. The predicted and observed rates of SS dileptons were compared as a function of kinematic properties and found to agree within 30%. The data driven estimate was also compared to a direct prediction from simulation and a similar level of agreement was reached.

The final category is a subdominant background in all SRs and corresponds to events where the charge of a lepton is incorrectly measured. Charge misidentification primarily occurs when an electron undergoes bremsstrahlung in the tracker material or in the beam pipe. Similarly to the tight-to-loose method, the number of SS lepton pairs where one of the leptons has its charge misidentified can be determined using the number of OS pairs and the knowledge of the charge misidentification rate. We use simulation to parameterize this rate as a function of \(p_{{\mathrm {T}}}\) and \(\eta \) for electrons and find values varying between \(10^{-5}\) (central electrons with \(p_{{\mathrm {T}}} \approx 20\,{\text {GeV}} \)) and \(5\times 10^{-3}\) (forward electrons with \(p_{{\mathrm {T}}} \approx 200\,{\text {GeV}} \)). To calibrate the charge misidentification rate, we exploit the fact that charge misidentification only has a small effect on the electron energy measurement in the calorimeter. As a result, electron pairs from \({\text {Z}}\) boson decays yield a sharp peak near the \({\text {Z}}\) mass even when one of the electrons has a misidentified charge. The SS dielectron invariant mass distributions in data and MC can then be used to derive a correction factor to the MC charge misidentification rate. Good agreement between data and MC is found in 2016, while the charge misidentification rate in simulation corresponding to 2017 and 2018 data needs to be scaled up by a factor of 1.4. Muon charge misidentification arises from a relatively large uncertainty in the transverse momentum at high momentum or from a poor quality track. The various criteria applied in this analysis on the quality of the muon reconstruction lead to a misidentification rate at least one order of magnitude smaller than for electrons according to simulation. The muon charge misassignment has also been studied using cosmic ray muons with \(p_{{\mathrm {T}}}\) up to several hundred \(\,{\text {GeV}}\), confirming the predictions from simulation [72]. It is therefore neglected. Correction factors are however applied to the simulation to account for a possible difference in the selection efficiency related to these criteria.

6 Systematic uncertainties

The predicted yields of signal and background processes are affected by several sources of uncertainty, summarized in Table 8. Depending on their source, they are treated as fully correlated or uncorrelated between the three years of data taking. Signal and background contributions estimated from simulation are affected by experimental uncertainties in the efficiency of the trigger, lepton reconstruction and identification [64, 73], the efficiency of \({\text {b}}\) tagging [69], the jet energy scale [67], the integrated luminosity [74,75,76]. An uncertainty is also assigned to the value of the inelastic cross section, which affects the pileup rate [77] and that can impact the description of the jet multiplicity or the \(p_{{\mathrm {T}}} ^{\text {miss}}\) resolution. Simulation is also affected by theoretical uncertainties, which are evaluated by varying the factorization and renormalization scales up and down by a factor of two, and by using different PDFs within the NNPDF3.0 and 3.1 PDF sets [35, 36, 78]. These uncertainties can affect both the overall yield (normalization) and the relative population (shape) across the SRs. Background normalization uncertainties are increased to 30%, either to account for the additional hadronic activity required (for \({\text {W}}{\text {Z}}\) and \({\text {W}}^{\pm } {\text {W}}^{\pm }\)) or to take into consideration recent measurements (for \(\hbox {t}{\bar{\hbox {t}}} {\text {W}}\), \(\hbox {t}{\bar{\hbox {t}}} {\text {Z}}\)) [79, 80]. The Rare and X\({\upgamma }\)backgrounds, which are less well understood experimentally and theoretically, are assigned a 50% uncertainty.

To account for possible mismodeling of the flavor of additional jets, an additional 70% uncertainty is applied to \(\hbox {t}{\bar{\hbox {t}}} {\text {W}}\), \(\hbox {t}{\bar{\hbox {t}}} {\text {Z}}\), and \(\hbox {t}{\bar{\hbox {t}}} {\text {H}} \) events produced in association with a pair of \({\text {b}}\) jets, reflecting the measured ratio of \(\hbox {t}{\bar{\hbox {t}}} {\text {b}}{\overline{\text {b}}} / \hbox {t}{\bar{\hbox {t}}} \text {jj}\) cross sections reported in Ref. [81].

As discussed in Sect. 5, the nonprompt lepton and charge misidentification backgrounds are estimated from CRs. The associated uncertainties include the statistical uncertainties in the CR yields, as well as the systematic uncertainties in the extrapolations from the CRs to the SRs, as described below. In the case of the nonprompt lepton background, we include a 30% uncertainty from studies of the closure of the method in simulation. Furthermore, the uncertainty in the measurement of the tight-to-loose ratio, because of the prompt lepton contamination, results in a 1–30% additional uncertainty in the background yields. The charge misidentification background is assigned a 20% uncertainty based on a comparison of the kinematic properties of simulated and data events in the \({\text {Z}}\rightarrow {\text {e}}^{+} {{\text {e}}}^{-} \) CR with one electron or positron having a misidentified charge.

In general, the systematic uncertainties with the largest impact on the expected limits defined below are related to the lepton identification and isolation scale factors, the cross section of the rare processes, and the \({\text {WZ}}\) background normalization.

Table 8 Summary of the sources of systematic uncertainty and their effect on the yields of different processes in the SRs. The first two groups list experimental and theoretical uncertainties assigned to processes estimated using simulation, while the last group lists uncertainties assigned to processes whose yield is estimated from the data. The uncertainties in the first group also apply to signal samples. Reported values are representative for the most relevant signal regions

7 Results and interpretation

The distributions of the variables used to define the SRs after the event selection are shown in Fig. 4. Background yields shown as stacked histograms in Figs. 45, and 6 are those determined following the prescriptions detailed in Sect. 5. The overall data yields exceed expectation by an amount close to the systematic uncertainty. However, no particular trend that is not covered by the uncertainties discussed in the previous sections, is seen in the distributions. The significance of the excess is of similar magnitude in all categories, with a maximum of around 2 standard deviations (s.d.) in the off-\({\text {Z}}\) ML category.

The results of the search, broken down by SR, are presented in Figs. 5 and 6, and are summarized in Table 9. No significant deviation with respect to the SM background prediction is observed. The largest excess of events found by fitting the data with the background-only hypothesis is in HH SR54, corresponding to a local significance of 2.6 s.d. Its neighboring bin, HH SR55, which is adjacent along the \(H_{\text {T}}\) dimension, has a deficit of events in the data corresponding to a significance of 1.8 s.d.

These results are then interpreted as experimental constraints on the cross sections for the signal models discussed in Sect. 2. For each model, event yields in all SRs are used to obtain exclusion limits on the production cross section at 95% confidence level (\({\text {CL}}\)) with an asymptotic formulation of the modified frequentist \(\text {CL}_\text {s}\) criterion [82,83,84,85], where uncertainties are incorporated as nuisance parameters and profiled [84]. This procedure takes advantage of the differences in the distribution of events amongst the SR between the various SM backgrounds and the signal considered. The normalizations of the various backgrounds are in particular allowed to float within their uncertainties in the global fit, resulting in several backgrounds (nonprompt lepton, \(\hbox {t}{\bar{\hbox {t}}} {\text {W}}/{\text {Z}}/{\text {H}} \) and rare processes) being pulled up by around 1 s.d. for most of the signal points considered, which are often characterized by a distinctive distribution of events across the SRs. This observation is consistent with the current measurements of \(\hbox {t}{\bar{\hbox {t}}} {\text {W}}\) and \(\hbox {t}{\bar{\hbox {t}}} {\text {Z}}\) processes performed by the ATLAS and CMS Collaborations [79, 80]. The limits obtained are then used together with the theoretical cross section calculations to exclude regions of SUSY parameter space.

Fig. 4
figure 4

Distributions of the main analysis variables after the event selection: \(H_{\text {T}}\), \(p_{{\mathrm {T}}} ^{\text {miss}}\), \(m_{\text {T}} ^{\text {min}}\), \(N_\text {jets}\), \(N_{{\text {b}}}\), and the charge of the SS pair, where the last bin includes the overflow (where applicable). The hatched area represents the total statistical and systematic uncertainty in the background prediction. The lower panels show the ratio of the observed event yield to the background prediction. The prediction for the SUSY model \(\mathrm {T}1{{\text {tttt}}}\) with \(m_{{\tilde{\text {g}}}}=1600\,{\text {GeV}} \) and \(m_{{\tilde{\chi }}^{0}_{1}}=600\,{\text {GeV}} \) is overlaid

Fig. 5
figure 5

Expected and observed SR yields for the HH, HL, LL signal categories. The hatched area represents the total statistical and systematic uncertainty in the background prediction

Fig. 6
figure 6

Expected and observed SR yields for the LM, on-\({\text {Z}}\) ML, off-\({\text {Z}}\) ML signal categories. The hatched area represents the total statistical and systematic uncertainty in the background prediction

Table 9 Expected background event yields, total uncertainties, and observed event yields in the SRs used in this search

Figure 7 shows observed and expected exclusion limits for simplified models of gluino pair production with each gluino decaying to off- or on-shell third-generation squarks. These models were introduced in Sect. 2 and denoted as \(\mathrm {T}1{{\text {tttt}}}\), \(\mathrm {T}5{{\text {ttbbWW}}}\), \(\mathrm {T}5{{\text {tttt}}}\), and \(\mathrm {T}5{{\text {ttcc}}}\). Similarly, Figs. 8 and 9 show the corresponding limits for \(\mathrm {T}5{{\text {qqqqWZ}}}\) and \(\mathrm {T}5{{\text {qqqqWW}}}\), with two different assumptions on the chargino mass. Note that the \(\mathrm {T}5{{\text {qqqqWZ}}}\) model assumes equal probabilities for the decay of the gluino into \({\tilde{\chi }}^{+}_{1}\), \({\tilde{\chi }}^{-}_{1}\), and \({\tilde{\chi }} _{2}^{0}\). The exclusion limits for \(\mathrm {T}6{{\text {ttWW}}}\) and \(\mathrm {T}6{{\text {ttHZ}}}\) are displayed in Figs. 10 and 11, respectively. In the \(\mathrm {T}6{{\text {ttHZ}}}\) model, the heavier top squark decays into a lighter top squark and a \({\text {Z}}\) or \({\text {H}}\) boson. The three sets of exclusion limits shown in Fig. 11 correspond to the branching fraction \({\mathcal {B}}({\tilde{\text {t}}} _{2}\rightarrow {\tilde{\text {t}}} _{1}{\text {Z}})\) having values of 0, 50, and 100%.

Finally, Fig. 12 shows observed and expected limits on the cross section of gluino pair production as a function of the gluino masses in the two RPV models described in Sect. 2. The observed and expected exclusions on the gluino mass are similar and reach 2.1 and \(1.7\,{\text {TeV}} \) for the \(\mathrm {T}1{{\text {qqqqL}}} \) and \(\mathrm {T}1{{\text {tbs}}}\) models, respectively.

The analysis sensitivity for the various models studied in Figs. 7, 8, 9, 10 and 11 is often driven by the event yields in a few SRs (off-\(\text {Z}\) ML21, HH53 and HH52), where a slight excess of data is observed. This in particular applies to the uncompressed mass regime, resulting in an observed limit weaker than the expected one by one or two s.d. In the compressed mass regime, however, other SRs can become dominant, for example when the hadronic activity becomes limited. This happens in the \(\mathrm {T}5{{\text {qqqqWZ}}}\) and \(\mathrm {T}5{{\text {qqqqWW}}}\) models where the gluino and the lightest neutralino present a limited mass splitting (the region close to the diagonal in the left plots of Figs. 8, 9). In those scenarios the on-\({\text {Z}}\) ML4 and HH3 SRs provide the best sensitivity, respectively. Additionally, if the intermediate chargino is nearly degenerate in mass with the lightest neutralino, both leptons become soft and LL SRs such as LL2 become relevant. Such a situation is encountered in the phase space region close to the diagonal in the right plots of Figs. 8 and  9. On-\({\text {Z}}\) SRs (especially on-\({\text {Z}}\) ML23) become important for models where an on-shell \({\text {Z}}\) boson is produced (bottom plot in Fig. 11). The limits on the RPV models presented in Fig. 12 are mostly driven by another set of SRs (HH62 and LM11, the latter becoming more relevant for lower masses).

Compared to the previous versions of the analysis [24, 25], the limits for the RPC models extend the gluino and squark mass observed and expected exclusions by up to \(200\,{\text {GeV}} \) because of the increase in the integrated luminosity and the corresponding re-optimization of SR definitions. These results also complement searches for gluino pair production conducted by CMS in final states with 0 or 1 lepton [86,87,88]. For the \(\mathrm {T}1{{\text {tttt}}}\) scenario, the expected sensitivity of this analysis suffers from a lower branching fraction that makes it uncompetitive in the uncompressed mass regime. However, for a nearly degenerate mass spectrum, the SM background becomes of higher importance and the presence of an SS lepton pair significantly reduces it, leading to a similar sensitivity. The constraints on the two RPV models that were not previously included demonstrate the sensitivity of the analysis to RPV scenarios. The final state is particularly well suited to study the \(\mathrm {T}1{{\text {qqqqL}}} \) model since no leptonic branching fraction penalty applies, resulting in exclusion limits on the gluino mass beyond 2.1 TeV, comparable to other results in fully hadronic final states [87, 88]. The limits obtained on the \(\mathrm {T}1{{\text {tbs}}}\) model are stronger than those previously obtained in the one-lepton channel based on the analysis of the 2016 dataset [89]. They are expected to remain competitive after an update with the full Run 2 dataset.

Model-independent limits are also set on the product of cross section, branching fraction, detector acceptance, and reconstruction efficiency, for the production of an SS lepton pair with at least two extra jets and \(H_{\text {T}} >300\,{\text {GeV}} \). For this purpose, we select events from the HH and LM categories and calculate limits as a function of minimum \(p_{{\mathrm {T}}} ^{\text {miss}}\) or \(H_{\text {T}}\) requirements starting at 300 and \(1400\,{\text {GeV}} \), respectively. In order to remove the overlap between the two conditions, events selected for the \(H_{\text {T}}\) scan must also satisfy \(p_{{\mathrm {T}}} ^{\text {miss}} <300\,{\text {GeV}} \). The corresponding limits are presented in Fig. 13.

Finally, in order to facilitate reinterpretations of our results, we present in Table 10 the expected and observed yields for a number of inclusive SRs. This procedure focuses on events with large \(H_{\text {T}}\), \(p_{{\mathrm {T}}} ^{\text {miss}}\), \(N_{{\text {b}}}\), and/or \(N_\text {jets}\), and the SRs are defined such that they typically lead to 5 to 10 expected background events. The last column in the table indicates the upper limit at 95% \({\text {CL}}\) on the number of BSM events in each SR.

Fig. 7
figure 7

Exclusion regions at 95% \({\text {CL}}\) in the \(m_{{\tilde{\chi }}^{0}_{1}}\) versus \(m_{{\tilde{\text {g}}}}\) plane for the \(\mathrm {T}1{{\text {tttt}}}\)  (upper left) and \(\mathrm {T}5{{\text {ttbbWW}}}\)  (upper right) models, with off-shell third-generation squarks, and the \(\mathrm {T}5{{\text {tttt}}}\)  (lower left) and \(\mathrm {T}5{{\text {ttcc}}}\) (lower right) models, with on-shell third-generation squarks. For the \(\mathrm {T}5{{\text {ttbbWW}}}\) model, \(m_{{\tilde{\chi }}^\pm _{1}} = m_{{\tilde{\chi }}^{0}_{1}} + 5\,{\text {GeV}} \), for the \(\mathrm {T}5{{\text {tttt}}}\) model, \(m_{{\tilde{\text {t}}}} - m_{{\tilde{\chi }}^{0}_{1}} = m_{{\text {t}}}\), and for the \(\mathrm {T}5{{\text {ttcc}}}\) model, \(m_{{\tilde{\text {t}}}} - m_{{\tilde{\chi }}^{0}_{1}} = 20\,{\text {GeV}} \) and the decay proceeds through \({\tilde{\text {t}}} \rightarrow {\text {c}}{\tilde{\chi }}^{0}_{1} \). The right-hand side color scale indicates the excluded cross section values for a given point in the SUSY particle mass plane. The solid black curves represent the observed exclusion limits assuming the approximate-NNLO+NNLL cross sections [46,47,48,49,50,51, 58] (thick line), or their variations of \(\pm 1\) standard deviations (s.d.) (thin lines). The dashed red curves show the expected limits with the corresponding \(\pm 1\) s.d. and \(\pm 2\) s.d. uncertainties. Excluded regions are to the left and below the limit curves

Fig. 8
figure 8

Exclusion regions at 95% \({\text {CL}}\) in the plane of \(m_{{\tilde{\chi }}^{0}_{1}}\) versus \(m_{{\tilde{\text {g}}}}\) for the \(\mathrm {T}5{{\text {qqqqWZ}}}\) model with \(m_{{\tilde{\chi }}^\pm _{1}}=0.5(m_{{\tilde{\text {g}}}} + m_{{\tilde{\chi }}^{0}_{1}})\) (left) and with \(m_{{\tilde{\chi }}^\pm _{1}} = m_{{\tilde{\chi }}^{0}_{1}} + 20\,{\text {GeV}} \) (right). The notations are as in Fig. 7

Fig. 9
figure 9

Exclusion regions at 95% \({\text {CL}}\) in the plane of \(m_{{\tilde{\chi }}^{0}_{1}}\) versus \(m_{{\tilde{\text {g}}}}\) for the \(\mathrm {T}5{{\text {qqqqWW}}}\) model with \(m_{{\tilde{\chi }}^\pm _{1}}=0.5(m_{{\tilde{\text {g}}}} + m_{{\tilde{\chi }}^{0}_{1}})\) (left) and with \(m_{{\tilde{\chi }}^\pm _{1}} = m_{{\tilde{\chi }}^{0}_{1}} + 20\,{\text {GeV}} \) (right). The notations are as in Fig. 7

Fig. 10
figure 10

Exclusion regions at 95% \({\text {CL}}\) in the plane of \(m_{{\tilde{\chi }}^\pm _{1}}\) versus \(m_{{\tilde{\text {b}}} _{1}}\) for the \(\mathrm {T}6{{\text {ttWW}}}\) model with \(m_{{\tilde{\chi }}^{0}_{1}}=50\,{\text {GeV}} \). The notations are as in Fig. 7

Fig. 11
figure 11

Exclusion regions at 95% \({\text {CL}}\) in the plane of \(m({\tilde{\text {t}}} _{1})\) versus \(m({\tilde{\text {t}}} _{2})\) for the \(\mathrm {T}6{{\text {ttHZ}}}\) model with \(m({\tilde{\text {t}}} _{1})-m({\tilde{\chi }}^{0}_{1})=175\,{\text {GeV}} \). The three exclusions represent \({\mathcal {B}}({\tilde{\text {t}}} _{2}\rightarrow {\tilde{\text {t}}} _{1}{\text {Z}})\) of 0, 50, and 100%, respectively. The notations are as in Fig. 7

Fig. 12
figure 12

Upper limits at 95% \({\text {CL}}\) on the cross section for RPV gluino pair production with each gluino decaying into four quarks and one lepton (\(\mathrm {T}1{{\text {qqqqL}}} \), left), and each gluino decaying into a top, bottom, and strange quarks (\(\mathrm {T}1{{\text {tbs}}}\), right)

Fig. 13
figure 13

Upper limits at 95% \({\text {CL}}\) on the product of cross section, detector acceptance, and selection efficiency, \(\sigma \! {\mathcal {A}} \epsilon \), for the production of an SS lepton pair with at least two jets, as a function of the minimum \(p_{{\mathrm {T}}} ^{\text {miss}}\) threshold, when \(H_{\text {T}} >300\,{\text {GeV}} \) (left), or the minimum \(H_{\text {T}}\) threshold, when \(p_{{\mathrm {T}}} ^{\text {miss}} <300\,{\text {GeV}} \) (right)

Table 10 Inclusive SR definitions, expected background yields and uncertainties, and observed yields, as well as the observed 95% \({\text {CL}}\) upper limits on the number of BSM events contributing to each region. No uncertainty in the signal acceptance is assumed in calculating these limits. A dash (—) indicates that a particular selection is not required

8 Summary

A sample of events with two same-sign or at least three charged leptons (electrons or muons) produced in association with several jets in proton-proton collisions at \(13\,{\text {TeV}} \), corresponding to an integrated luminosity of \(137{\,{\text {fb}}^{-1}} \), has been studied to search for manifestations of physics beyond the standard model. The data are found to be consistent with the standard model expectations. The results are interpreted as limits on cross sections at 95% confidence level for the production of new particles in simplified supersymmetric models, considering both R parity conserving and violating scenarios. Using calculations for these cross sections as functions of particle masses, the limits are translated into lower mass limits that are as large as \(2.1\,{\text {TeV}} \) for gluinos and \(0.9\,{\text {TeV}} \) for top and bottom squarks, depending on the details of the model. The results extend the gluino and squark mass observed and expected exclusions by up to \(200\,{\text {GeV}} \), compared to the previous versions of this analysis. Finally, to facilitate further interpretations of the search, model-independent limits are provided as a function of the missing transverse momentum and the scalar sum of jet transverse momenta in an event, together with the background prediction and data yields in a set of simplified signal regions.