Search for new phenomena in three-or four-lepton events in pp collisions at

A search with minimal model dependence for physics beyond the Standard Model in events featuring three or four charged leptons (3 (cid:3) and 4 (cid:3) , (cid:3) = e , μ ) 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 √ s = 13 TeV and recorded with the ATLAS detector, corresponding to the full Run 2 dataset of 139 fb − 1 . The 3 (cid:3) and 4 (cid:3) 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 signiﬁcant 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 (pp) collisions with exactly three or four leptons (3 and 4 , where = e, μ in this paper).An example is supersymmetry (SUSY), where neutralino and chargino production [10] yields three leptons and a neutrino through an intermediate W Z 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 E-mail address: atlas .publications@cern .ch.
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 AT-LAS 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 crosssection, 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 pp dataset collected by the ATLAS experiment during the 2015-2018 data-taking period is used, cor-responding 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 Z -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 W Z and Z Z 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 pp collisions at √ s = 8 TeV [29].
Compared to that search, the current analysis uses a larger dataset collected at √ s = 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] 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 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 z-axis along the beam pipe.The x-axis points from the IP to the centre of the LHC ring, and the y-axis points upwards.Cylindrical coordinates (r, φ) are used in the transverse plane, φ being the azimuthal angle around the z-axis.The pseudorapidity is defined in terms of the polar angle θ as η = − ln tan(θ/2).Angular distance is measured in units of R ≡ ( η)2 + ( φ) 2 .
that reduces the accepted event rate to 1 kHz on average depending on the data-taking conditions.

Data and simulated events
This analysis uses pp collision data at a centre-of-mass energy of √ s = 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 (ee, eμ and μμ) are used.
The transverse momentum (p 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, they are 18 and 8 GeV in 2015 and 22 and 8 GeV in 2016-18. 2 For the mixed-flavour trigger, the p 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 p 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. 3The dominant SM backgrounds are the production of two vector bosons decaying leptonically: W Z for the 3 final states and Z Z for the 4 final states.Subleading prompt-lepton backgrounds that contribute are triboson production, and processes which include at least one top quark: t t X ( X = W , Z , H ), t Z , t tW W and t tt t.
All diboson and triboson (V V and V V V , where V = W , Z ) production, including off-shell production, was simulated with the Sherpa 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 Sherpa authors.The matrix element calculations were matched and merged with the Sherpa 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 V V jj ( j = jet) events were generated at LO.This contribution includes Higgs boson production through vector-boson fusion (with H → Z Z ).It also includes triboson processes where one boson decays hadronically, including V H → V V V → V V jj.Loop-induced production of Z Z 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 gg → H → Z Z .
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 OpenLoops library [48][49][50].This process includes only on-shell fully leptonic decays.
The production of t tV , t tW W , t Zq and t tt t events was modelled using the MadGraph5_aMC@NLO v2.3 [51] (v2.2 for t tW W ) generator, while t t H events [52] were modelled using the PowhegBox [53][54][55][56] v2 generator with the h damp parameter set to 1.5 m top [57].Events were generated at NLO for t tV and t t H.The t Zq, t tW W and t tt t 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 PYTHIA8.210 (8.230 for t Zq) [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 (H ±± ) particles decaying leptonically.The simplified Type-III seesaw model was included in the Mad-Graph5_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 H ±± events were generated at LO using the left-right-symmetry package of PYTHIA8.186, which provides the H ±± scenario described in Ref. [28].Implementation details are described in Ref. [21].The parton shower was provided by PYTHIA8.230(PYTHIA8.186)for the Type-III seesaw (H ±± ) model using the NNPDF2.3loPDF set and the A14 tune [59].Cross-sections for both samples were normalised to NLO.
The ATLAS detector simulation [61] employing the Geant4 [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 pp events onto hard-scatter events.These were generated with PYTHIA8 [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 p 2 T of its associated tracks, which must include at least two with p 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 p T > 25 GeV and |η| < 2.47.For the purposes of matching tracks to the primary vertex, the track impact parameters4 d 0 and z 0 must satisfy |d 0 |/σ (d 0 ) < 5 (3) for electrons (muons) and |z 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 p T > 300 GeV, the High-p T working point [66] is used, which places a tighter requirement on the number of muon spectrometer hits, ensuring optimal momentum resolution for highly energetic muons.
Both the electron and muon candidates are required to be isolated in the ID.To determine isolation, a cone is placed around the object's track, with an opening angle which is the smaller of R = 0.2 (0.3) for electrons (muons) and 10 GeV/p T, .The scalar p T sum of all tracks (excluding the lepton itself) within this cone, I R , must satisfy I R /p T, < 0.06.For muon candidates with p T > 50 GeV, the opening angle of the isolation cone is always R = 0.2.Electron candidates must also pass a calorimeter-based isolation requirement of I R /p T, < 0.06, this time taking the sum of calorimeter energy deposits as I R , within a cone of R = 0.2.
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-k t algorithm [69,70] with a radius parameter R = 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 p T > 20 GeV and |η| < 2.5.The jet vertex tagger (JVT) [72] is used to test jets that have p 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 p T is rejected.Jets are rejected if they are within R = 0.2 (for overlap removal, the pseudorapidity in R is substituted with the rapidity, defined as y = − ln E+p z E−p z ) 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 R = 0.4 of any remaining jets are removed.
The missing transverse momentum (E 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 p 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 E miss T and the presence of an on-Z lepton pair: a same-flavour and oppositely charged (SFOC) lepton pair with a dilepton mass within 10 GeV of the Z -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

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 Z -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 Z -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.

Region
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 E miss where the E 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-Z lepton pair, and are called on-Z and off-Z SRs respectively.No selection based on charge and flavour is made for these SRs beyond the SFOC pair needed for the on-Z region.Certain heavy BSM particles are expected to decay into lepton pairs without first decaying to an intermediate Z -boson (e.g.Ref. [28]).Off-Z SRs are expected to be sensitive to such signals while excluding the main prompt-lepton background contribution (leptonically decaying W Z and Z Z vector-boson pairs, which are likely to be on-Z ).Furthermore, on-Z SRs include the few events where three leptons combine into multiple valid on-Z pairs.SRs are further split according to the distribution of the invariant mass (m inv ) of all leptons in the event.Four divisions are established to construct the 3 regions: 0-200 GeV, 200-400 GeV, 400-600 GeV, and >600 GeV.Two divisions are established for the 4 regions: 0-400 GeV and >400 GeV.This leads to 22 SRs in total.
Two CRs are defined: a 3 CR for the W Z background and a 4 CR for the Z Z background.The W Z control region requires an on-Z lepton pair and a third off-Z lepton which has a transverse mass (m T , defined as 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 Z + jets and t t processes 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 E miss T of the event (E miss T < 25 GeV for electrons, E 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 ( j) with p T > 35 GeV and φ(μ, j) > 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 p 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 eeμ + eμμ events are designated as on-Z and off-Z 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-Z and off-Z VRs consist of different ratios of fake-lepton sources: the on-Z VR is more sensitive to Z +jets events than the off-Z region, while the reverse is true for t t events.Furthermore, both VRs have a substantial W Z contribution which can be validated.In the off-Z VR the proper modelling of the W Z MC simulation of off-shell Z decays is confirmed.For this VR, a SFOC pair of off-Z leptons is required.Both validation regions target, through a m T requirement of m T ( , E 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 W -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 pileup 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 E miss T computation.Further uncertainties affecting the E 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-Z region, which is insensitive to E miss T ) and 4% for the 4 , on-Z , 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 W +jets and Z +jets contributions).The fake-lepton background estimate is compared with an estimate obtained when using fake factors that consider events with and without b-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 E miss T upper bound in the fake-lepton estimation regions, and by a p 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 m inv bin and region.
For the t t X ( X = W , Z , H ) 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 t tW and t t Z processes (Sherpa [39] for t tW , MadGraph5_aMC@NLO [51] interfaced with Herwig 7 [85] for t t Z ) 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 t Zq, t tW W and t tt t, 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 E miss T and m 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 m 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 (W Z for the 3 SRs, Z Z 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-Z SRs with m inv < 400 GeV); for the two 3 SRs The ATLAS Collaboration Physics Letters B 824 (2022) 136832 Fig. 1.Comparison between data and prediction for the m inv distribution of the W Z (left) and Z Z (right) control regions after the fit to the data.The rightmost bin is inclusive and contains all events with m inv > 500 GeV.'Fakes' refers to the fake-lepton background.The hatched grey area shows the combination of all uncertainties in the analysis.

Fig. 2. Comparison between data and prediction for the E miss
T distribution of the W Z (left) and Z Z (right) control regions after the fit to the data.The rightmost bin is inclusive and contains all events with E miss T > 300 GeV.'Fakes' refers to the fake-lepton background.The hatched grey area shows the combination of all uncertainties in the analysis.
where E miss T < 50 GeV and m inv > 400 GeV; and for the 3 , off-Z , E miss T > 50 GeV SR with m 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: N S .The same W Z and Z Z 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 N 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(n|μ) for the single SR and multiple control region bins.The predicted number of events in each region is expressed in terms of the W Z and Z Z background normalisation factors k W Z and k Z Z , the nuisance parameters and, in the case of the SR, N S .For each systematic uncertainty, indexed by l, the likelihood formula is multiplied by a constraint term, which is a standard Gaussian response function G(0|θ l , σ l ) with mean and standard deviation equal to the central value of the nuisance parameter and its uncertainty, respectively.The likelihood therefore takes the form: Comparison between data and prediction for the m inv distribution of the on-Z (left) and off-Z (right) validation regions after the fit to the data.The rightmost bin is inclusive and contains all events with m inv > 500 GeV.'Fakes' refers to the fake-lepton background.The hatched grey area shows the combination of all uncertainties in the analysis.

Fig. 4.
Comparison between data and prediction in each signal region of this analysis after the profile likelihood fit has been performed.'Fakes' refers to the fake-lepton background.The hatched grey area shows the combination of all uncertainties in the analysis.
where the index j runs over all background contributions other than W Z and Z Z .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).

Likelihood fit to data
Figs. [1][2][3] show comparisons between data and predictions for CRs and VRs after performing a likelihood fit to the CRs.The values of the normalisation factors are obtained from a binned CR-only fit (i.e. its likelihood formula does not include a term for any signal region) to the distributions in Fig. 1.The normalisation factors are found to be 0.98 ± 0.07 for the W Z background, consistent with the generator cross-section, and 1.05 ± 0.09 for the Z Z background, which is consistent with previous measurements of the event yield of on-shell Z Z decays [75].Furthermore, the VRs show good agreement between data and background, indeed validating the off-shell W Z modelling and the understanding of the fakelepton background.
Fig. 4 shows all SRs with the observed data and post-fit background yields.The event yields obtained by comparing the measured data with the expected background are shown, after performing the fitting procedure, in Table 2 for the 3 signal regions, and in Table 3 for the 4 signal regions.
The assumed number of signal events N 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

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.backgrounds further; the post-fit background values for every SR are therefore always the same as those of the CR-only fit, with N 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 ( NS ) and, from the parabolic behaviour of the log-likelihood around its maximum, its associated uncertainty ( NS ).From these values a signif- icance can be computed, defined as Z = NS / NS so that negative significance is associated with negative signal yields.The signifi-cances for all SRs are given in Table 4.In this analysis, all SRs have a significance below |Z | < 3.

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] 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.

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 A × ). 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].

Model
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 (N 95 ) by the integrated luminosity L of the collected data: σ vis = N 95 /L.These limits are given for all SRs in Fig. 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 modelindependent 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 H ±± 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 m L = 400 GeV and 7.5 +3.1 −1.8 fb for m L = 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 m H ±± = 300 GeV and 0.14 +0.13 −0.07 fb for m H ±± = 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 .

Conclusion
In this paper, a model-independent search targeting final states with three or four light leptons is presented using 139 fb −1 of √ s = 13 TeV pp 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 Z -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 crosssections.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].

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.and NSF, United States of America.In addition, individual groups and members have received support from BCKDF, CANARIE, Compute Canada and CRC, Canada; COST, ERC, ERDF, Horizon 2020 and Marie Skłodowska-Curie Actions, European Union; Investissements d'Avenir Labex, Investissements d'Avenir Idex and ANR, France; DFG and AvH Foundation, Germany; Herakleitos, Thales and Aristeia programmes co-financed by EU-ESF and the Greek NSRF, Greece; BSF-NSF and GIF, Israel; Norwegian Financial Mechanism 2014-2021, Norway; NCN and NAWA, Poland; La Caixa Banking Foundation, CERCA Programme Generalitat de Catalunya and PROMETEO and GenT Programmes Generalitat Valenciana, Spain; Göran Gustafssons Stiftelse, Sweden; The Royal Society and Leverhulme Trust, United Kingdom.

T.
This motivates a selection of SRs separated by a E miss T cut.A threshold of E miss T = 50 GeV was chosen, which splits the phase space into regions where the E miss T originates mostly from detector resolution effects and regions of 40 < m T < 80 GeV, which captures leptons originating from a Wboson decay.The Z Z control region requires four leptons to form two on-Z lepton pairs.The SRs are separated from the W Z CR through their flavour composition or by requiring m T > 80 GeV, and from the Z Z CR by vetoing events with 4 on-Z 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.

Fig. 3 .
Fig. 3.Comparison between data and prediction for the m inv distribution of the on-Z (left) and off-Z (right) validation regions after the fit to the data.The rightmost bin is inclusive and contains all events with m inv > 500 GeV.'Fakes' refers to the fake-lepton background.The hatched grey area shows the combination of all uncertainties in the analysis.

1 − 1 . 3 Fig. 5 .
Fig.5.The expected and observed 95% CL upper limits for each signal region of this analysis expressed in terms of the visible cross-section σ vis .

Table 4
Local significance of the value of the parameter of interest for each signal region after performing the combined likelihood fit, defined as Z = NS / NS .
h i j k l n Also at Department of Physics, California State University, Fresno; United States of America.o Also at Department of Physics, California State University, Sacramento; United States of America.p Also at Department of Physics, King's College London, London; United Kingdom.q Also at Department of Physics, St. Petersburg State Polytechnical University, St. Petersburg; Russia.r Also at Department of Physics, University of Fribourg, Fribourg; Switzerland.s Also at Faculty of Physics, M.V. Lomonosov Moscow State University, Moscow; Russia.t Also at Faculty of Physics, Sofia University, 'St.Kliment Ohridski', Sofia; Bulgaria.u Also at Giresun University, Faculty of Engineering, Giresun; Turkey.v Also at Graduate School of Science, Osaka University, Osaka; Japan.w Also at Hellenic Open University, Patras; Greece.