Search for new phenomena in three- or four-lepton events in $pp$ collisions at $\sqrt{s} = 13$ TeV with the ATLAS detector

A search with minimal model dependence for physics beyond the Standard Model in events featuring three or four leptons ($3\ell$ and $4\ell$, $\ell = e,\mu$) is presented. The analysis aims to be sensitive to a wide range of potential new-physics theories simultaneously. This analysis uses data from $pp$ collisions delivered by the Large Hadron Collider at a centre-of-mass energy of $\sqrt{s} = 13$ TeV and recorded with the ATLAS detector, corresponding to the full Run 2 dataset of 139 fb$^{-1}$. The $3\ell$ and $4\ell$ phase space is divided into 22 event categories according to the number of leptons in the event, the missing transverse momentum, the invariant mass of the leptons, and the presence of leptons originating from a $Z$-boson candidate. These event categories are analysed independently for the presence of deviations from the Standard Model. No statistically significant deviations from the Standard Model predictions are observed. Upper limits for all signal regions are reported in terms of the visible cross-section.


Introduction
Despite the success of the Standard Model (SM) [1][2][3][4] in describing the interactions of elementary particles, there remain observations that suggest the existence of additional phenomena [5][6][7][8][9]. Many theories of physics beyond the Standard Model (BSM theories) have been proposed that feature final states in high-energy proton-proton ( ) collisions with exactly three or four leptons (3ℓ and 4ℓ, where ℓ = , in this paper). An example is supersymmetry (SUSY), where neutralino and chargino production [10] yields three leptons and a neutrino through an intermediate state, and where di-Higgs production [11] yields four or more leptons. Furthermore, enhanced flavour-changing decay cross-section of top quarks at the loop level may lead to anomalous production of 3ℓ final states with respect to the SM [12]. Multiple types of seesaw models can produce multilepton final states alongside neutrinos [13][14][15][16][17]. Adding an additional Higgs triplet to the SM Lagrangian [18,19] potentially leads to a doubly charged Higgs particle [20,21] which can decay into two leptons, leading to a four-lepton final state if produced in pairs. Theories predicting such a particle include left-right symmetric models [22,23], scalar singlet dark matter [24] and the Zee-Babu model [25]. Conclusive evidence for any of these BSM theories has thus far been elusive. Many dedicated analyses within the LHC experimental collaborations are being performed to search for such evidence. However, the vast number of theories means that it would be difficult to perform a dedicated analysis for each model. This motivates the desire to establish instead a search that does not rely on a specific model for its signal description (henceforth called 'model-independent'), which can cover a wide range of signatures to seek indicators of exotic physics.
The analysis presented here is committed to investigating a large phase space while making few prior assumptions about the nature of new-physics processes. As such, it is expected to be sensitive to a large number of signals that could be populating the ATLAS data, albeit with a lower sensitivity than a dedicated analysis could achieve. Partial overlap with dedicated analyses is expected, but since these are typically tuned to specific models they do not consider the full phase space. For instance, Ref. [26] also studies 4ℓ final states, but does not consider events with low four-lepton invariant masses and missing transverse momenta. The present analysis does not attempt to exploit very distinctive features, such as e.g. resonances in invariant-mass distributions, and instead entails a substantially more inclusive selection than is typical of searches for specific BSM theories, Nevertheless, a comparison with a few benchmark models indicates that its sensitivity is not greatly reduced compared to that of more dedicated searches.
The analysis aims to uncover evidence of BSM physics. Failing that, it can provide a set of upper limits on the visible cross-section, which can be reinterpreted as upper limits on BSM models of interest. The upper limits on two benchmark models derived in this fashion, a Type-III seesaw model [27] and a doubly-charged Higgs model [28] are compared with those obtained using dedicated analyses.
For this analysis, the full dataset collected by the ATLAS experiment during the 2015-2018 data-taking period is used, corresponding to an integrated luminosity of 139 fb −1 delivered by the LHC. Events featuring exactly three or four charged leptons are categorised into signal regions based on the invariant mass of the leptons, the missing transverse momentum and the presence of a lepton pair compatible with originating from a -boson decay. The observations in such regions are individually used to probe for the presence of a BSM signal. Control regions are established to extract a normalisation of the most prominent SM backgrounds, which are leptonically decaying and diboson pairs, and for the estimation of the contribution from interactions producing lepton candidates from heavy-flavor hadrons decays or hadronic particles misidentified as leptons. In each region, a SM-only hypothesis is compared with a hypothesis assuming the SM plus an additional number of BSM events as a free parameter.
A previous general multilepton search with the ATLAS detector was performed using 20.3 fb −1 of collisions at √ = 8 TeV [29]. Compared to that search, the current analysis uses a larger dataset collected at √ = 13 TeV and assigns 4ℓ events to multiple separate regions, but does not consider hadronic -lepton decays. A strategy for a general search was outlined in Ref.
[30] using 3.2 fb −1 of data at √ = 13 TeV. That search uses fewer data events and a coarser background estimation, but offers a broader selection of final states, including multilepton final states, and performs tests with additional variables. Furthermore, a similar search, also testing multilepton final states, has been performed by the CMS Collaboration with 137 fb −1 of collisions at √ = 13 TeV [31].

ATLAS detector
The ATLAS experiment [32] at the LHC is a multipurpose particle detector with a forward-backward symmetric cylindrical geometry and a near 4 coverage in solid angle. 1 It consists of an inner tracking detector (ID) surrounded by a thin superconducting solenoid providing a 2 T axial magnetic field, electromagnetic and hadron calorimeters, and a muon spectrometer. The inner tracking detector covers the pseudorapidity range | | < 2.5. It consists of silicon pixel (with the insertable B-layer installed before Run 2 [33,34]), silicon microstrip, and transition radiation tracking detectors. In the range | | < 3.2, hermetic lead/liquid-argon (LAr) sampling calorimeters provide electromagnetic (EM) energy measurements with high granularity. The central region, | | < 1.8, is additionally instrumented with a thin LAr presampling detector to correct for energy losses in the inactive material in front of the detector. A steel/scintillator-tile hadronic calorimeter covers the central pseudorapidity range | | < 1.7. The endcap and forward regions are instrumented with LAr calorimeters for both the EM and hadronic energy measurements up to | | = 4.9. The muon spectrometer surrounds the calorimeters. It consists of three large superconducting air-core toroidal magnets with eight coils each. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. The muon spectrometer includes three layers of precision tracking chambers, allowing precise muon momentum measurements up to | | = 2.7, and fast detectors for triggering up to | | = 2.4. A two-level trigger system [35] is used to select events. The first-level trigger is implemented in hardware and uses a subset of the detector information to accept events at a rate below 100 kHz. This is followed by a software-based trigger that reduces the accepted event rate to 1 kHz on average depending on the data-taking conditions.

Data and simulated events
This analysis uses collision data at a centre-of-mass energy of √ = 13 TeV, collected by the ATLAS detector and corresponding to a total integrated luminosity of 139 fb −1 . Only data recorded during stable beam conditions with all ATLAS detector subsystems operational [36] have been included. Dilepton triggers [37,38] covering all lepton flavour combinations ( , and ) are used. The transverse momentum ( T ) requirements of these triggers depend on the data-taking period. For the dielectron trigger, these requirements are 12 GeV in 2015, 17 GeV in 2016, and 24 GeV in 2017-18. For the dimuon trigger, 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the centre of the detector and the -axis along the beam pipe. The -axis points from the IP to the centre of the LHC ring, and the -axis points upwards. Cylindrical coordinates ( , ) are used in the transverse plane, being the azimuthal angle around the -axis. The pseudorapidity is defined in terms of the polar angle as = − ln tan( /2). Angular distance is measured in units of they are 18 and 8 GeV in 2015 and 22 and 8 GeV in 2016-18. 2 For the mixed-flavour trigger, the T  requirement is 17 GeV for the electron and 14 GeV for the muon. These trigger choices correspond to the  dilepton triggers with the lowest T requirements available during each data-taking year. Expected event rates due to SM processes that can result in 3ℓ and 4ℓ final states were estimated using a combination of Monte Carlo (MC) event generation and data-driven techniques. Event generators based on MC methods were used to estimate the total expected contributions from SM processes producing only prompt leptons. 3 The dominant SM backgrounds are the production of two vector bosons decaying leptonically: for the 3ℓ final states and for the 4ℓ final states. Subleading prompt-lepton backgrounds that contribute are triboson production, and processes which include at least one top quark: ( = , , ), ,¯and¯¯.
All diboson and triboson ( and , where = , ) production, including off-shell production, was simulated with the S 2.2.2 [39] generator. The NNPDF3.0nnlo set of PDFs was used [40], along with the dedicated set of tuned parton-shower parameters developed by the S authors. The matrix element calculations were matched and merged with the S parton shower [41] based on Catani-Seymour dipole factorisation [42,43] using the MEPS@NLO prescription [44][45][46][47]. Diboson events were generated at next-to-leading-order (NLO) accuracy in QCD for up to one additional parton and at leading-order (LO) accuracy for two and three additional parton emissions. Electroweak ( = jet) events were generated at LO. This contribution includes Higgs boson production through vector-boson fusion (with → ). It also includes triboson processes where one boson decays hadronically, including → → . Loop-induced production of events via gluon-gluon fusion was simulated using matrix elements accurate at LO for up to one additional parton emission for both the fully leptonic and semileptonic final states. This contribution includes → → .
Triboson events were generated at NLO for the inclusive process and at LO for up to two additional parton emissions. The virtual QCD corrections were provided by the O L library [48][49][50]. This process includes only on-shell fully leptonic decays.
The production of¯,¯, and¯¯events was modelled using the M G 5_aMC@NLO v2.3 [51] (v2.2 for¯) generator, while¯events [52] were modelled using the P B [53][54][55][56] v2 generator with the ℎ damp parameter set to 1.5 top [57]. Events were generated at NLO for¯and¯. The , and¯¯processes were modelled at LO with their cross-sections normalised to NLO predictions [51]. The NNPDF3.0nlo [40] PDF was used. The events were interfaced to P 8.210 (8.230 for ) [58] using a set of tuned parameters called the A14 tune [59] and the NNPDF2.3lo [40] PDF set. Collectively, these events are referred to as the top-quark background.
Additional processes involving Higgs boson production, apart from those mentioned above, have been found to yield negligible contributions and are not considered explicitly For model-specific interpretation of the analysis, signal samples were generated for Type-III seesaw model heavy leptons and doubly charged Higgs ( ±± ) particles decaying leptonically. The simplified Type-III seesaw model was included in the M G 5_aMC@NLO v2.3.3 generator at LO with an implementation using FeynRules [60], described in Ref. [27]. Implementation details are described in Ref. [17]. Samples of ±± events were generated at LO using the left-right-symmetry package of P 8.186, which provides the ±± scenario described in Ref. [28]. Implementation details are described in Ref. [21]. The parton shower was provided by P 8.230 (P 8.186) for the Type-III seesaw ( ±± ) model using the NNPDF2.3lo PDF set and the A14 tune [59]. Cross-sections for both samples were normalised to NLO.
The ATLAS detector simulation [61] employing the G 4 [62] framework was used to model the detector response in MC events. The effect of pile-up was incorporated into the simulation by overlaying additional inelastic events onto hard-scatter events. These were generated with P 8 [58] using the A3 tune [63] and the MSTW2008LO [64] PDF set. Events were simulated with discrete values for the expected mean number of interactions and then weighted to match the distribution that is observed per bunch crossing in data.

Object selection
All events in the analysis are required to have a primary vertex, defined as the vertex with the highest value of 2 T of its associated tracks, which must include at least two with T > 0.5 GeV. Electron and muon candidates are required to originate from the primary vertex.
Requirements common to both electron and muon candidates are T > 25 GeV and | | < 2.47. For the purposes of matching tracks to the primary vertex, the track impact parameters 4 0 and 0 must satisfy | 0 |/ ( 0 ) < 5 (3) for electrons (muons) and | 0 sin | < 0.5 mm. Furthermore, electron and muon candidates are subjected to identification criteria, for which multiple working points are provided [65,66]. The identification is performed using quality cuts where each working point offers a different trade-off between the rate of false positives and false negatives delivered by the algorithm.
Electron candidates are reconstructed using energy clusters measured in the EM calorimeter matched to reconstructed tracks [65]. The identification for the nominal selection in this analysis is based on a combination of detailed tracking and calorimeter information combined into a likelihood discriminant. The Tight [65] working point is used. The range 1.37 < | | < 1.52 has a significant amount of non-sensitive material in front of the calorimeter, and is therefore excluded.
Muon candidates are reconstructed by combining measurements in the ID and the muon spectrometer [67]. For this analysis, the Medium identification working point is used for most muon candidates [66], which requires an ID track matched with multiple muon spectrometer precision hits. For muon candidates with The constituents for jet reconstruction are identified by combining measurements from both the ID and the calorimeter using a particle-flow algorithm [68]. Jet candidates are reconstructed from these particle-flow objects using the anti-algorithm [69,70] with a radius parameter = 0.4. The jet energy scale (JES) and resolution (JER) [71] are corrected to particle level using MC simulation. Jets are furthermore required to have T > 20 GeV and | | < 2.5. The jet vertex tagger (JVT) [72] is used to test jets that have T < 60 GeV and | | < 2.4 to suppress those originating from pile-up.
Objects found to have very collinear tracks are considered to be overlapping. Overlaps are resolved through a sequence of rules. This procedure prevents double-counting of particles interacting with different parts of the detector, and provides an optimal classification of these particles. If a muon candidate is found to have a shared ID track with an electron candidate, the electron candidate is rejected. If two electron candidates have shared ID tracks, the one with the lower T is rejected. Jets are rejected if they are within Δ = 0.2 (for overlap removal, the pseudorapidity in Δ is substituted with the rapidity, defined as = − ln + − ) of a lepton candidate, except if the candidate is a muon and three or more collinear tracks are found. Subsequently, lepton candidates that are within Δ = 0.4 of any remaining jets are removed.
The missing transverse momentum ( miss T ) [73] in a given reconstructed event is computed as a combination of a hard term, the magnitude of the negative vector sum of the T of all reconstructed leptons and jets, and a soft term, computed from the momenta of inner-detector tracks that are not matched to any of the selected objects but do originate from the primary vertex.

Analysis strategy
Selected events are separated into different categories, referred to as regions, to maximise the sensitivity to a relatively broad range of potential new phenomena. Signal regions (SRs) are defined as regions to be probed for the presence of such signatures. Criteria that separate these SRs are the number of leptons, the miss T and the presence of an on-lepton pair: a same-flavour and oppositely charged (SFOC) lepton pair with a dilepton mass within 10 GeV of the -boson mass of 91.2 GeV. 5 Control regions (CRs) are defined so as to be dominated by particular SM processes which have been well-studied in the ATLAS experiment [74,75]. The full list of regions with their selection criteria, apart from the splitting into different invariant-mass ranges, is shown in Table 1. All regions are orthogonal to each other, so no event is assigned to more than one region. The CRs are used to extract the normalisation of the main background processes from data and to constrain the size of the systematic uncertainties of the analysis. Validation regions (VRs) are used to confirm that the predictions for SM background processes are well-modelled.
A large group of BSM models predict the existence of at least one additional heavy lepton beyond the SM (e.g. Ref. [27]), either charged or neutral. Such theories often feature final states with one or more neutrinos, which due to being invisible to the detector translates into a non-zero miss T . This motivates a selection of SRs separated by a miss T cut. A threshold of miss T = 50 GeV was chosen, which splits the phase space into regions where the miss T originates mostly from detector resolution effects and regions where the miss T is likely to be due to objects invisible to the detector.
SRs are also categorised according to the presence or absence of at least one on-lepton pair, and are called on-and off-SRs respectively. No selection based on charge and flavour is made for these SRs beyond the SFOC pair needed for the on-region. Certain heavy BSM particles are expected to decay into lepton pairs without first decaying to an intermediate -boson (e.g. Ref. [28]). Off-SRs are expected to Table 1: Overview of the regions defined for this analysis. The -symbols indicate that no requirements are made on the variable for that particular region. Additional requirements on SRs, described in the Other column, veto events that are used in the CR/VRs from entering into the SRs. The -pairs column denotes the number of non-overlapping lepton pairs that are same-flavour and oppositely charged and have a dilepton invariant mass within 10 GeV of the -boson mass of 91.2 GeV. The off-flavour ℓ is the lepton in the 3ℓ event that has a different flavour from the other two leptons; cuts requiring an off-flavour ℓ are not applied if all three leptons are of the same flavour.

Validation regions
3ℓ On- be sensitive to such signals while excluding the main prompt-lepton background contribution (leptonically decaying and vector-boson pairs, which are likely to be on-). Furthermore, on-SRs include the few events where three leptons combine into multiple valid on-pairs. Two CRs are defined: a 3ℓ CR for the background and a 4ℓ CR for the background. The control region requires an on-lepton pair and a third off-lepton which has a transverse mass ( T , defined as GeV, which captures leptons originating from a -boson decay. The control region requires four leptons to form two on-lepton pairs. The SRs are separated from the CR through their flavour composition or by requiring T > 80 GeV, and from the CR by vetoing events with 4 on-leptons. Both CRs are used to constrain the two normalisation factors of their corresponding backgrounds. These normalisation factors are free parameters in the statistical analysis. A data-driven technique is used to estimate backgrounds with at least one fake lepton, referred to as the fake-lepton background, in the SRs, CRs and VRs. Fake leptons are either non-prompt leptons or hadrons misidentified as leptons by the detector. The primary sources of such events are the + jets andp rocesses which have two prompt leptons and at least one fake lepton. The yield of fake-lepton background events is measured separately for electrons and muons using the fake-factor method, which is described in Ref. [29]. Dedicated regions containing a single lepton candidate are established using data collected by the single-lepton triggers. Selection requirements for these regions are imposed to ensure a large number of events with fake leptons, in order to reduce the statistical uncertainty of this contribution. Requirements are based on the miss T of the event ( miss T < 25 GeV for electrons, miss T < 40 GeV for muons) and on the number of jets in the event (≥ 1 for the electron region, ≥ 2 for the muon region). For the muon, there must also be at least one jet ( ) with T > 35 GeV and Δ ( , ) > 2.7, called the tag jet. For each region, an adjacent 'anti-ID' region is established, with orthogonal selection criteria of the identification and isolation algorithms. The anti-ID regions admit non-isolated electrons and muons, as well as isolated but Loose [65] electrons, but veto events which satisfy the nominal selection criteria. The anti-ID selection is tuned for a large number of fake leptons, while the stringent working points of the nominal selection suppress these fake leptons, to ensure the analysis is robust against fake-lepton contamination. A ratio is computed from the event rates of these two regions. This ratio is called the fake factor and is parameterised as a function of the T and of the lepton. An anti-ID region is also established for each signal, control and validation region of the analysis. These anti-ID regions use the same selection criteria, except that one or more of their leptons passes the alternative identification and isolation requirements. Using the fake factor, the yield of fake-lepton background events is extrapolated from each anti-ID region to its nominal counterpart. Prompt backgrounds, estimated through MC methods and normalised to their theoretical cross sections, are subtracted from observed data in the fake-factor estimation regions prior to calculation of the fake factor, and in the anti-ID regions prior to extrapolation.
Certain subselections of + events are designated as on-and off-VRs. These are used to check that the computed fake factors transfer correctly from the regions where they are calculated to the regions in which they are applied. The on-and off-VRs consist of different ratios of fake-lepton sources: the on-VR is more sensitive to +jets events than the off-region, while the reverse is true for¯events. Furthermore, both VRs have a substantial contribution which can be validated. In the off-VR the proper modelling of the MC simulation of off-shell decays is confirmed. For this VR, a SFOC pair of off-leptons is required. Both validation regions target, through a T requirement of T (ℓ, miss T ) < 40 GeV, a third lepton that is likely to be fake. Only mixed-flavour final states are selected for these VRs so that the choice of the third lepton, assumed to be the fake lepton (or the lepton due to -boson decay), is unambiguous.

Systematic uncertainties
Systematic uncertainties affect the precision of the predicted background contributions. Two classes of systematic uncertainties are defined: detector-related uncertainties, referred to as 'experimental', and uncertainties in MC modelling of the processes, referred to as 'theoretical'.

Experimental uncertainties
Multiple experimental uncertainties have been considered for this analysis, although only a small number of them have a significant impact on the results. These uncertainties are discussed below.
The uncertainty in the combined 2015-2018 integrated luminosity [76] is 1.7%, obtained using the LUCID-2 detector [77] for the primary luminosity measurements. Uncertainties in the reweighting procedure applied to the simulation to bring its pile-up multiplicity distribution into agreement with that in the data are also included, ranging between 0.5% and 1%.
For the leptons, uncertainties due to the measured momentum resolution and scale are taken into account [65,78]. Uncertainties in the reconstruction, identification and isolation efficiency scale factors that are used to correct for the difference between the MC simulation and data are also included. The impact of this set of uncertainties on the expected yield for the signal regions varies between 1% and 2%.
Sets of uncertainties in the jet energy scale and resolution are also included. These were derived from information taken from test-beam data, LHC collision data and simulation [71,79]. These uncertainties are small for all signal regions, with an impact between 0.2% and 0.9%.
Uncertainties associated with the above objects are propagated to an uncertainty in the hard term of the miss T computation. Further uncertainties affecting the miss T that are included are uncertainties in the offset and resolution of the soft term [80]. The impact of these ranges between no impact (for the 4ℓ, off-region, which is insensitive to miss T ) and 4% for the 4ℓ, on-, miss T > 50 GeV region.
Several systematic uncertainties on the fake factors are considered. First, there is an uncertainty due to limited number of events in the single-lepton region where the fake factors are calculated. Then, there is an uncertainty in the MC modelling of the dominant background contributions in the single-lepton region (the +jets and +jets contributions). The fake-lepton background estimate is compared with an estimate obtained when using fake factors that consider events with and without -jets [81], with the difference between the two estimates taken as an uncertainty. Finally, two uncertainties are included to address the bias caused by imposing a miss T upper bound in the fake-lepton estimation regions, and by a T requirement on the tag jet in the fake-muon estimation regions. These uncertainties are estimated by varying the requirements on these variables upwards and downwards by 10 GeV. The impact of the fake-factor uncertainties on the total background prediction ranges between 0.1% (for the 4ℓ regions, where there are very few fake-lepton events) and 1.6%.

Theoretical uncertainties
Theoretical uncertainties affect the MC-based background estimate of the multiboson and top-quark backgrounds. The main theoretical uncertainties considered for this analysis originate from the missing higher orders in the perturbative expansion of the partonic cross-section, from PDF uncertainties and the choice of PDF, and from the uncertainty in the strong coupling constant ( s ). The analysis follows the PDF4LHC recommendations [82] for the computation of these uncertainties. These uncertainties are uncorrelated between different background contributions. Other uncertainties such as matching and merging uncertainties, hadronisation and parton-shower uncertainties are not included as this analysis is not directly sensitive to jets.
For the diboson and triboson processes, the contribution of missing higher-order diagrams is estimated by observing the differences in the cross-section prediction when varying the renormalisation scale r and factorisation scale f . These scales are independently varied upwards and downwards by a factor of two [83], leaving out predictions where the terms are scaled in opposite directions. This leads to a total of seven scale variations. The total uncertainty is taken as the envelope of all variations, picking the variation with the largest value in each individual inv bin and region.
For the¯( = , , ) contributions to the top-quark background, the uncertainty due to missing higher orders is estimated in the same way as for the diboson and triboson uncertainties, using the envelope of the seven variations of r and f , while the PDF uncertainty is taken as the standard deviation of 100 replica variations. The impact of uncertainties in s is taken from Ref. [84]. Predictions from alternative generators for¯and¯processes (S [39] for¯, M G 5_aMC@NLO [51] interfaced with H 7 [85] for¯) are used in assessing an uncertainty due to the choice of generator; this uncertainty is found to have no impact on the final result. For the rare top-quark processes ,¯and¯¯, the scale and PDF uncertainties are taken as those associated with the computed NLO cross-section values reported in Ref. [51]. A more precise estimation of these uncertainties is considered unnecessary due to the minute contribution of these rare top-quark processes to the total background yield of this analysis.
The impact of the scale uncertainties on the total background estimate ranges between 5% and 15%; the SRs with high miss T and inv requirements are at the high end of this range as the impact of higher-order diagrams is especially large there. The impact of the PDF uncertainty is around 2%-3%, with higher values (up to 6%) in higher inv bins. The impact of the s uncertainty is 1%-2%.
In general, the dominant sources of systematic uncertainty for this analysis are the theoretical uncertainties. Of these, the r and f scale uncertainty of the diboson backgrounds ( for the 3ℓ SRs, of the 4ℓ SRs) has the largest effect. However, many SRs are still statistically limited. This is the case for the 4ℓ SRs (except for the 4ℓ, on-SRs with inv < 400 GeV); for the two 3ℓ SRs where miss T < 50 GeV and inv > 400 GeV; and for the 3ℓ, off-, miss T > 50 GeV SR with inv > 600 GeV.

Statistical analysis and results
Each of the 22 SRs is treated in an individual, together with the CRs, cut-and-count experiment, and a statistical analysis is performed independently for each region. For each of these analyses, the parameter of interest is the number of signal events in the corresponding SR: S . The same and control regions are used for all statistical analyses as well.
This analysis employs a maximum-likelihood technique, using the profile likelihood ratio (see, e.g., Ref. [86]) to estimate S while also accounting for the various systematic uncertainties affecting the background predictions, which enter the likelihood expression as nuisance parameters (NP) ì . The likelihood for each SR is the product of Poisson probability terms P ( | ) for the single SR and multiple control region bins. The predicted number of events in each region is expressed in terms of the and background normalisation factors and , the nuisance parameters and, in the case of the SR, S . For each systematic uncertainty, indexed by , the likelihood formula is multiplied by a constraint term, which is a standard Gaussian response function G(0| , ) with mean and standard deviation equal to the central value of the nuisance parameter and its uncertainty, respectively. The likelihood therefore takes the form: where the index runs over all background contributions other than and . Many sources of uncertainty affect the predictions in both the SR and the CRs. The correlations between predictions in the different regions are accounted for by a common dependence on the associated nuisance parameter(s). The impact of these may vary between control and signal regions, particularly for the dominant systematic uncertainties (which are the diboson scale uncertainties).   Table 2 for the 3ℓ signal regions, and in Table 3 for the 4ℓ signal regions.

Likelihood fit to data
The assumed number of signal events S is allowed to float freely. It is determined independently for each SR and is not required to be positive. Given that the SR is not binned, it does not add further degrees of freedom and hence does not constrain the backgrounds further; the post-fit background values for every SR are therefore always the same as those of the CR-only fit, with S adjusted to make up for the difference between the total background and the observed data in the SR.
For each SR a fit finds the number of signal events (ˆS) and, from the parabolic behaviour of the log-likelihood around its maximum, its associated uncertainty (ΔˆS). From these values a significance can be computed, defined as =ˆS/ΔˆS so that negative significance is associated with negative signal yields. The significances for all SRs are given in Table 4. In this analysis, all SRs have a significance below | | < 3.         Table 2: Summary of the event yields for all background contributions to the 3ℓ signal regions after the combined likelihood fit has been performed. The observed data for each signal region are also given. 'Fakes' refers to the fake-lepton background. The 'Total' row gives the sum of event counts of all individual Standard Model contributions. The 'Signal' row gives the remaining discrepancy between the sum of backgrounds and observed data. The total uncertainty in the event yield is given for each background contribution, for the total Standard Model prediction and for the best-fit value of the signal contribution. It should be noted that the uncertainty in the total background prediction cannot be obtained as a sum in quadrature of the uncertainties in its contributions, due to statistical correlations.   Table 3: Summary of the event yields for all background contributions to the 4ℓ signal regions after the combined likelihood fit has been performed. The observed data for each signal region are also given. 'Fakes' refers to the fake-lepton background. The 'Total' row gives the sum of event counts of all individual Standard Model contributions. The 'Signal' row gives the remaining discrepancy between the sum of backgrounds and observed data. The total uncertainty in the event yield is given for each background contribution, for the total Standard Model prediction and for the best-fit value of the signal contribution. It should be noted that the uncertainty in the total background prediction cannot be obtained as a sum in quadrature of the uncertainties in its contributions, due to statistical correlations.

Visible cross-section limits
No significant excess was found in any of the signal regions. Limits on the number of signal events are set.
The CL S method [87] is used to ascertain upper limits in the signal regions. Assumptions made about the test statistic are based on the works of Wilks [88] and Wald [89]. Specifically, it is assumed that the test statistic asymptotically approaches a 2 distribution with one degree of freedom for a large number of events [90].
Expected and observed 95% CL upper limits are presented for the visible cross-section vis , which is calculated by dividing the upper limit on the total number of events ( 95 ) by the integrated luminosity L of the collected data: vis = 95 /L. These limits are given for all SRs in Figure 5. Visible cross-section limits in this figure can be reinterpreted as limits on specific physics models as long as the selection efficiency and acceptance of the model (including any uncertainties in these values) for a specific SR definition used in this analysis is known. By dividing the visible cross-section limits given here by this efficiency and acceptance, upper limits on the cross-section can be derived from this analysis.

Model-specific limits
The results are interpreted for two particular BSM models previously studied by dedicated analyses. The purpose of this interpretation is to compare the results obtained by this model-independent search with those from a dedicated search. Studied models are the Type-III seesaw model described in Ref.
[17] and the search for doubly charged Higgs boson production described in Ref. [21]. The chosen parameters of the models studied in this section are at parity with the referenced analyses, but only two representative particle mass hypotheses are chosen for each models: 400 GeV and 700 GeV for the heavy lepton in the Type-III seesaw model, and 300 GeV and 500 GeV for the ±± particles. These masses are chosen to cover a range of masses corresponding to the simulated models available and the published limits available for comparison.
The signal regions and background predictions remain the same as described earlier in this Letter. Using these models, all of the signal regions in these analyses are studied to find the region with the best limit-setting power, judged by the value of the expected limit. The best signal region and the corresponding expected and observed limits found by this analysis are given in Table 5. To convert the limits from the model-independent analysis to limits on the cross-section of the signal model considered, only a correction for acceptance effects and selection efficiency is applied. This procedure ignores uncertainties in these quantities; however, the experimental uncertainties on the expected signal yields in the most sensitive regions are not larger than a few percent.
The Type-III seesaw model analysis [17] presents an expected 95% CL cross-section upper limit of 22 +8.5 −6.4 fb for = 400 GeV and 7.5 +3.1 −1.8 fb for = 700 GeV for the full Run 2 dataset with an integrated luminosity of 139 fb −1 , although it only tests dilepton final states. These limits are more stringent than those derived with the analysis presented in this paper that correspond to 41 and 12 fb, respectively. The doubly-charged Higgs boson analysis [21] presents an expected 95% CL cross-section upper limit for the four-lepton final state of 0.16 +0.14 −0.07 fb for ±± = 300 GeV and 0.14 +0.13 −0.07 fb for ±± = 500 GeV. These limits are comparable to those when using the best limit of a single SR of this analysis. However, the upper limits in Ref. [21] are obtained using only 2015 and 2016 ATLAS data, corresponding to an integrated luminosity of 36.1 fb −1 . Table 5: Expected and observed cross-section exclusion limits at 95% CL for representative mass values of the two selected models. Also the most sensitive bin, which was used to obtain these limits for each case, is listed, along with the signal acceptance times efficiency in this region (denoted by × ).

Model
Mass

Conclusion
In this paper, a model-independent search targeting final states with three or four light leptons is presented using 139 fb −1 of √ = 13 TeV collision data collected by the ATLAS detector at the LHC from 2015 to 2018. The analysis offers a wide coverage of the 3ℓ and 4ℓ phase space. The measured data of 3ℓ and 4ℓ events is tested for potential indicators of physics beyond the Standard Model. By categorising the targeted phase space according to the number of leptons, the missing transverse momentum, the presence of a lepton pair originating from a -boson decay, and the invariant mass of the leptons in the event, a total of 22 signal regions are defined. Each signal region is analysed independently using a profile likelihood fit. Control regions were established to extract the normalisations of the primary background processes.
No significant deviations from the Standard Model expectation are found in the data. In the absence of a detected signal, upper limits at the 95% CL are provided in terms of the visible cross-sections. The expected upper limits can be interpreted by dedicated analyses as long as the efficiency and acceptance of their signal model in a particular signal region is known. The analysis is interpreted using simulated signal models for heavy leptons from the Type-III seesaw mechanism [17] and a doubly charged Higgs boson model [21].
(Taiwan), RAL (UK) and BNL (USA), the Tier-2 facilities worldwide and large non-WLCG resource providers. Major contributors of computing resources are listed in Ref. [91].
[10] ATLAS Collaboration, Search for chargino-neutralino production with mass splittings near the electroweak scale in three-lepton final states in