Search for physics beyond the standard model in events with jets and two same-sign or at least three charged leptons in proton-proton collisions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=13\,{\text {TeV}} $$\end{document}s=13TeV

A data sample of events from proton-proton collisions with at least two jets, and two isolated same-sign or three or more charged leptons, is studied in a search for signatures of new physics phenomena. The data correspond to an integrated luminosity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$137{\,{\text {fb}}^{-1}} $$\end{document}137fb-1 at a center-of-mass energy of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$13\,{\text {TeV}} $$\end{document}13TeV, collected in 2016–2018 by the CMS experiment at the LHC. The search is performed using a total of 168 signal regions defined using several kinematic variables. The properties of the events are found to be consistent with the expectations from standard model processes. Exclusion limits at 95% confidence level are set on cross sections for the pair production of gluinos or squarks for various decay scenarios in the context of supersymmetric models conserving or violating R parity. The observed lower mass limits are as large as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.1\,{\text {TeV}} $$\end{document}2.1TeV for gluinos and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.9\,{\text {TeV}} $$\end{document}0.9TeV for top and bottom squarks. To facilitate reinterpretations, model-independent limits are provided in a set of simplified signal regions.


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 (pp) 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 pairproduced gluinos or squarks can contain multiple W or Z bosons, with the potential to have at least one pair of SS W bosons. Such a decay chain is realized, for example in e-mail: cms-publication-committee-chair@cern.ch gluino pair production, when a gluino decays into a top quarkantiquark pair and a neutralino, or into a pair of quarks and a chargino that subsequently decays into a 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 pp collisions at √ s = 13 TeV recorded by the CMS detector at the CERN LHC in 2016-2018, corresponding to an integrated luminosity of 137 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 τ 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 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.

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 ttW, ttZ and WZ contributions) is simulated with the MadGraph5_amc@nlo 2.2.2 (2.4.2) [27][28][29] generator for 2016 (2017 and 2018) conditions. An exception is the WZ 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 W boson pairs and other very rare SM processes, are generated at leading order (LO) accuracy with Mad-Graph5_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 W ± W ± , 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 MadGraph5_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 ISR J ). The reweighting factors are extracted from a study of the light-flavor jet multiplicity in dilepton tt events. They vary between 0.92 and 0.77 for N ISR J 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 (  and where several leptons can arise in the final state. Production cross sections are calculated at approximate next-tonext-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 b quarks and up to four W bosons are shown in Fig. 1. In these models, the gluino decays to the lightest squark (g →qq), which in turn decays to same-flavor (q → qχ 0 1 ) or different-flavor (q → q χ ± 1 ) quarks. The chargino (χ ± 1 ) decays to a W boson and a neutralino (χ 0 1 ) viaχ ± 1 → W ±χ 0 1 , where theχ 0 1 is taken to be the lightest stable SUSY particle and escapes detection.
The first scenario, denoted by T1tttt and displayed in Fig. 1a, includes an off-shell top squark (t) leading to the three-body decay of the gluino,g → ttχ 0 1 , resulting in events with four W bosons and four b quarks. Figure 1b presents a similar model (T5ttbbWW) where the gluino decay results in a chargino that further decays into a neutralino and a W boson. The model shown in Fig. 1c (T5tttt) is the same as T1tttt except that the intermediate top squark is on-shell. The mass splitting between thet and theχ 0 1 is taken to be m˜t − mχ0 where m t is the top quark mass. This choice maximizes the kinematic differences between this model and T1tttt, and also corresponds to one of the most challenging regions of parameter space for the observation of thẽ t → tχ 0 1 decay since the neutralino is produced at rest in the top squark rest frame. The decay chain of Fig. 1d (T5ttcc) is identical to that of T5tttt except that thet decay involves a cquark. In Fig. 1e, the decay process includes a virtual lightflavor squark, leading to three-body decays ofg → qq χ ± 1 org → qq χ 0 2 , with a resulting signature of two W bosons, two Z bosons, or one of each (the case shown in Fig. 2e), and four light-flavor jets. This model, T5qqqqWZ, with a resulting signature of one W boson and one Z boson, is studied with two different assumptions for the chargino mass: mχ± 1 = 0.5(mg + mχ0 1 ), and mχ± 1 = mχ0 1 + 20 GeV, producing on-and off-shell bosons, respectively. The model is also considered with the assumption of decays to two W bosons exclusively (T5qqqqWW). Figure 2a shows a model of bottom squark production with subsequent decay ofb 1 → tχ ± 1 , yielding two b quarks and four W bosons. This model, T6ttWW, is considered as a function of the the lightest bottom squark,b 1 , andχ ± 1 masses. Theχ 0 1 mass is fixed to be 50 GeV, causing two of the W bosons to be produced off-shell when theχ ± 1 mass is less than approximately 130 GeV. Figure 2b displays a model similar to T6ttWW, but with top squark pair production and a subsequent decay oft 2 →t 1 H/Z, witht 1 → tχ 0 1 , producing signatures with two H bosons, two Z bosons, or one of each. In this model, T6ttHZ, theχ 0 1 mass is fixed such that m(t 1 ) − m(χ 0 1 ) = m t . The R parity violating decays considered in this analysis are T1qqqqL (Fig. 3a) and T1tbs (Fig. 3b). In T1qqqqL, the gluino decays to the lightest squark (g →qq), which in turn decays to a quark (q → qχ 0 1 ), but decays with theχ 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 T1tbs, each gluino decays into three different SM quarks (a top, a bottom, and a strange quark).

The CMS detector and event reconstruction
The central feature of the CMS detector is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant 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 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 kHz before data storage.
The reconstructed vertex with the largest value of summed physics-object squared-transverse-momentum is taken to be the primary 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 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 calorime-ters 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 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 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 b quarks are identified as btagged jets using a deep neural network algorithm, DeepCSV [69], with a working point chosen such that the efficiency to identify a b jet is 55-70% for a jet p T between 20 and 400 GeV. The misidentification rate for a light-flavor jet is 1-2% in the same jet p T range.
The vector p miss T 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 miss T . The scalar p T sum of all jets in an event is referred to as H T .

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 T > 40 GeV and two light SS leptons with p T > 15 GeV (10 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 T within a cone centered on the lepton direction and whose radius decreases with its p T , the ratio of the p T of the lepton to that of the closest jet, and the relative p 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 T > 25 GeV, increasing as a function of p T and reaching the maximum value for p T > 60 GeV. In the lowmomentum regime, 15 < p T < 25 GeV for electrons and 10 < p T < 25 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 T > 30 GeV, and it is around 92% for muons.
In order to reduce backgrounds from the decays of c-and b-hadrons or from the Drell-Yan process, we reject events with same-flavor lepton pairs with invariant mass (m ) less than 12 GeV, where leptons are reconstructed with a looser set of requirements compared to the nominal selection. Furthermore, events containing a lepton pair with m < 8 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 T of the leptons forming the SS pair: high-high if both have p T > 25 GeV, low-low if both have p T < 25 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 T thresholds on the leading (subleading) lepton in the 17-23 (8)(9)(10)(11)(12) GeV range, and dilepton triggers with no isolation requirements, a lower p T threshold of 8 GeV, an invariant mass condition m > 8 GeV to reject low mass resonances, and with a minimum H T in the range of 300−350 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 T requirements.
Six exclusive categories are then defined as follows: same as on-Z ML but without a 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 (T1tttt, T5ttbbWW, T5tttt, T1tttt, T5qqqqWW) with a large mass splitting between the gluino and the lightest neutralino. 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 Z bosons in the final state such as T5qqqqWZ and T6ttHZ 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 T1qqqqL where no genuine p miss T is expected. Various SRs are constructed based on the jet multiplicity N jets , the b-tagged jet multiplicity N b , H T , p miss T , the charge of the SS pair, and m min T , which is defined below. The m min T 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 p miss T , except for the on-Z ML category where we only consider the transverse mass computed using the leptons not forming the Z candidate. It is characterized by a kinematic cutoff for events where p miss T only arises from the leptonic decay of a single 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 WZ, ttW or SS WW. The SRs corresponding to each category, HH, HL, LL, LM, on-Z ML, and off-Z ML, are summarized in Tables 2, 3, 4, 5, 6 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 Table 2 The SR definitions for the HH category. Charge-split regions are indicated with (++) and (-- SR57 N jets = 5 or 6 SR58 N jets = 5 or 6 SR59 N jets = 5 or 6  Table 4 The SR definitions for the LL category. All SRs in this category require N jets ≥ 2. There are 8 regions in total. Quantities are specified in units of GeV where applicable. Table 5 The SR definitions for the LM category. All SRs in this category require p miss T < 50 GeV and H T > 300 GeV. The two high-H T regions are split only by N jets , resulting in 11 regions in total. Quantities are specified in units of GeV where applicable.
≥2 SR7 Table 6 The SR definitions for the on-ZML category. All SRs in these categories require N jets ≥ 2. Regions marked with † are split by m min T = 120 GeV, with the high-m min T region specified by the second SR label. There are 23 regions in total. Quantities are specified in units of GeV where applicable. Table 7 The SR definitions for the off-Zcategory. All SRs in these categories require N jets ≥ 2. Regions marked with † are split by m min T = 120 GeV, with the high-m min T region specified by the second SR label. There are 21 regions in total. Quantities are specified in units of GeV where applicable.
any SR has relative statistical uncertainties typically smaller than unity.

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 miss T or H T and in most of the ML SRs with a Z candidate. The main contributions arise from the production of a WZ or an SS Wpair, 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  or of a tt pair in association with a W, Z or H boson. The event yields for these processes are estimated individually. In contrast, the expected event yields from other rare processes (including ZZ, triple boson production, tWZ, tZq, tttt, and double parton scattering) are summed up into a single contribution denoted as "Rare". Processes including a genuine photon, such as Wγ, Zγ, ttγ, and tγ, are also considered and grouped together. They are referred to as "Xγ". 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 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 miss T or low m min T (except for the on-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 miss T , and low m min T . This sample is contaminated by prompt leptons from W boson decays. The contamination is estimated from the m min T 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 η and p corr T . The p corr T variable is defined as the sum of the lepton p 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 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 tt and 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 tightto-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 T and η for electrons and find values varying between 10 −5 (central electrons with p T ≈ 20 GeV) and 5 × 10 −3 (forward electrons with p T ≈ 200 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 Z boson decays yield a sharp peak near the 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 T up to several hundred 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.

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 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 miss T 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 WZ and W ± W ± ) or to take into consideration recent measurements (for ttW, ttZ) [79,80]. The Rare and Xγ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 ttW, ttZ, and ttH events produced in association with a pair of b jets, reflecting the measured ratio of ttbb/ttjj 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

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. 4, 5, 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-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 pre-  Fig. 7 diction 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 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 (CL) with an asymptotic formulation of the modified frequentist CL s criterion [82-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, ttW/Z/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 ttW and ttZ 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. 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 T1tttt, T5ttbbWW, T5tttt, and T5ttcc. Similarly, Figs. 8 and 9 show the corresponding limits for T5qqqqWZ and T5qqqqWW, with two different assumptions on the chargino mass. Note that the T5qqqqWZ model assumes equal probabilities for the decay of the gluino intoχ + 1 ,χ − 1 , andχ 0 2 . The exclusion limits for T6ttWW and T6ttHZ are displayed in Figs. 10 and 11, respectively. In the T6ttHZ model, the heavier top squark decays into a lighter top squark and a Z or H boson. The three sets of exclusion limits shown in Fig. 11 correspond to the branching fraction B(t 2 →t 1 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 TeV for the T1qqqqL and T1tbs 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-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 T5qqqqWZ and T5qqqqWW 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-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-Z SRs (especially on-Z ML23) become important for models where an on-shell 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 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-88]. For the T1tttt 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 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 T > 300 GeV. For this purpose, we select events from the HH and LM categories and calculate limits as a function of minimum p miss T or H T requirements starting at 300 and 1400 GeV, respectively. In order to remove the overlap between the two conditions, events selected for the H T scan must also satisfy p miss T < 300 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 T , p miss T , N b , and/or N 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% CL on the number of BSM events in each SR.

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 TeV, corresponding to an integrated luminosity of 137 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.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: Extended results
Tables 11, 12, 13, 14, 15 and 16, corresponding to Figs. 5 and 6, show background predictions per process within each signal region.     Appendix B: Top five SRs for several representative models Table 17 presents the top five SRs for several representative models, ranked based on the largest values of N sig. / √ N bkg. + N sig. , where N sig. and N bkg. are the signal and total background yields in each SR, respectively.