Search for direct stau production in events with two hadronic τ -leptons in ﬃﬃ s p = 13 TeV pp collisions with the ATLAS detector

A search for the direct production of the supersymmetric partners of τ -leptons (staus) in final states with two hadronically decaying τ -leptons is presented. The analysis uses a dataset of pp collisions corresponding to an integrated luminosity of 139 fb − 1 , recorded with the ATLAS detector at the Large Hadron Collider at a center-of-mass energy of 13 TeV. No significant deviation from the expected Standard Model background is observed. Limits are derived in scenarios of direct production of stau pairs with each stau decaying into the stable lightest neutralino and one τ -lepton in simplified models where the two stau mass eigenstates are degenerate. Stau masses from 120 GeV to 390 GeV are excluded at 95% confidence level for a massless lightest neutralino.

In SUSY models, the sector of sparticles with only electroweak interactions contains charginos (χ AE i , i ¼ 1, 2 in order of increasing masses), neutralinos (χ 0 j , j ¼ 1, 2, 3, 4 in order of increasing masses), charged sleptons (l), and sneutrinos (ν). Charginos and neutralinos are the mass eigenstates formed from linear superpositions of the superpartners of the charged and neutral Higgs bosons and electroweak gauge bosons. The sleptons are the superpartners of the charged leptons and are referred to as left or right (l L orl R ) depending on the chirality of their SM partners. The slepton mass eigenstates are a mixture ofl L andl R , and are labeled asl k (k ¼ 1, 2 in order of increasing mass). In this work, the scalar superpartner of the left-handed τ-lepton (the stau-left,τ L ) and right-handed τ-lepton (the stau-right,τ R ) are assumed to be mass degenerate.
Final states with τ-leptons originating from stau decays are of particular interest for SUSY searches. Models with light staus can lead to a dark-matter relic density consistent with cosmological observations [11] and light sleptons in general could play a role in the coannihilation of neutralinos [12,13]. Sleptons are expected to have masses of order 100 GeV in gauge-mediated [14][15][16] and anomalymediated [17,18] SUSY breaking models.
In some scenarios the direct production rate of sleptons can greatly exceed the production rate of squarks and gluinos at the Large Hadron Collider (LHC). In the simplified models studied in this paper, the lightest neutralino is the LSP and is purely the superpartner of the U(1) gauge field (the bino) and not admixed with the superpartner of the SU(2) gauge field (wino). Further, the two staus are assumed to be mass-degenerate. The stau-left and stau-right decay each with a 100% branching fraction into a binolike neutralino and a τ-lepton. All sparticles other than those explicitly mentioned here are assumed to be inaccessible at the LHC energy. This paper focuses on the direct production of a stau pair, leading to the final state illustrated in Fig. 1.
Signal events are characterized by the presence of exactly two τ-leptons and large missing transverse momentum due to the undetected neutrinos, as a result of the τ-lepton decays, and lightest neutralinos. Final states with exactly two hadronically decaying τ-leptons (τ → hadrons ν τ ) are considered.
The search described in this document uses a dataset of ffiffi ffi s p ¼ 13 TeV pp collisions collected with the ATLAS detector from 2015 to 2018 at the LHC. At LEP [19][20][21][22][23][24], searches set a lower limit of 86.6 GeV at 95% confidence level (C.L.) on the mass of promptly decaying staus. Similar searches were performed previously by the ATLAS and CMS Collaborations using data collected at ffiffi ffi s p ¼ 8 TeV in 2012 at the LHC [25,26], with the ATLAS result excluding a singleτ mass hypothesis of 109 GeV decaying to a τ-lepton and massless neutralino. CMS has performed a similar search using an integrated luminosity of 77.2 fb −1 of ffiffi ffi s p ¼ 13 TeV pp collisions and has excluded degenerate production ofτ pairs with mass up to 150 GeV for a nearly massless neutralino [27].

II. ATLAS DETECTOR
The ATLAS detector [28] is a multipurpose particle physics detector with forward-backward symmetric cylindrical geometry and nearly 4π coverage in solid angle. 1 It features an inner tracking detector (ID) surrounded by a 2 T superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer (MS). The ID covers the pseudorapidity region jηj < 2.5 and consists of a silicon pixel detector, a silicon microstrip detector, and a transition radiation tracker. One significant upgrade for the ffiffi ffi s p ¼ 13 TeV running period is the presence of the insertable B-layer [29,30], an additional pixel layer close to the interaction point which provides high-resolution hits at a small radius to improve the tracking and vertex reconstruction performance. The calorimeters are composed of high-granularity liquid-argon (LAr) electromagnetic calorimeters with lead, copper, or tungsten absorbers (in the pseudorapidity region jηj < 3.2) and a steelscintillator hadronic calorimeter (for jηj < 1.7). The end cap and forward regions are instrumented with LAr calorimeters for both the electromagnetic and hadronic measurements up to jηj ¼ 4.9. The MS surrounds the calorimeters and consists of three large superconducting air-core toroidal magnets, each with eight coils, a system of precision tracking chambers (jηj < 2.7), and detectors for triggering (jηj < 2.4). A two-level trigger system is used to select events for recording [31].

III. DATA AND SIMULATED EVENT SAMPLES
After the application of beam, detector, and data quality requirements, the analyzed dataset corresponds to an integrated luminosity of 139 fb −1 of pp collision data recorded from 2015 to 2018 at ffiffi ffi s p ¼ 13 TeV. The uncertainty in the combined 2015-2018 integrated luminosity is 1.7% [32], obtained using the LUCID-2 detector [33] for the primary luminosity measurements. The average number of interactions per bunch crossing (hμi) for this dataset ranges from about 10 up to 60, with a mean value of 34.
Monte Carlo (MC) simulated event samples were used to estimate the SUSY signal yields and to aid in evaluating the SM backgrounds. Generated SM events were passed through a detailed detector simulation [34] based on GEANT4 [35]. SUSY MC events were processed with a fast detector simulation that parameterizes the response of the electromagnetic and hadronic calorimeters but uses GEANT4 for the other detectors. All simulated events were overlaid with multiple pp collisions (pileup) simulated with the soft strong-interaction processes of PYTHIA 8.186 [36] using the A3 set of tuned parameters [37] and the NNPDF23LO [38] parton distribution function (PDF) set. The simulated events were reconstructed using the same algorithms as the data, and were reweighted so that the distribution of the expected number of collisions per bunch crossing matched the one in the data.

A. Simulated background samples
Events with Z=γ Ã → llðl ¼ e; μ; τÞ and W → lν produced with accompanying jets (including jets initiated by heavy flavor quarks) were generated with SHERPA 2.2.1 [39,40]. Matrix elements (ME) were calculated for up to two additional partons at NLO and four additional partons at leading order (LO), using the Comix [41] and OPENLOOPS [42,43] generators and merged with the SHERPA parton shower (PS) [44] using the ME þ PS@NLO prescription [40]. The NNPDF3.0NNLO [45] PDF set was used in conjunction with a dedicated PS tuning developed by the SHERPA authors. The W=Z þ jets events were normalized using their next-to-next-to-leading-order (NNLO) cross sections [46].
Fully leptonically and semileptonically decaying diboson samples (VV ¼ WW=WZ=ZZ) were simulated with the SHERPA 2.2.1 and 2.2.2 [39] generator at NLO. In this setup, multiple matrix elements were matched and merged with the SHERPA parton shower based on Catani-Seymour dipole factorization [44] using the MEPS@NLO prescription [40,[47][48][49]. The virtual QCD corrections for matrix elements at NLO accuracy were provided by the OPENLOOPS library [43]. Samples were generated using FIG. 1. A diagram illustrating the pair production of staus and subsequent decay into a two-τ-lepton final state with missing transverse momentum from the neutralinos. 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector, and the z-axis along the beam line. The x-axis points from the IP to the center 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. Observables labeled transverse refer to the projection into the x-y plane. The pseudorapidity is defined in terms of the polar angle θ by η ¼ − ln tanðθ=2Þ. the NNPDF3.0NNLO set, along with the dedicated set of tuned parton-shower parameters developed by the SHERPA authors.
The production of top-quark pairs (tt) and single top quarks in the Wt and s-channels was performed with POWHEG-BOX v2 [50][51][52][53], with the NNPDF2.3LO [54] PDF set at NLO in the ME calculations and the ATLAS underlying-event tune A14 [55]. Electroweak t-channel single-top-quark events were generated using the POWHEG-BOX v2 event generator. The PS, fragmentation, and the underlying event were simulated using PYTHIA 8.186 with the NNPDF2.3LO PDF set and a corresponding set of A14 tuned parameters. The top-quark mass was set to 172.5 GeV. The tt sample is normalised to the cross section prediction at next-to-next-to-leading order (NNLO) in QCD including the resummation of next-to-next-to-leading logarithmic (NNLL) soft-gluon terms calculated using TOP++2.0 [56][57][58][59][60][61][62]. The cross section for single-top-quark was computed for the Wt-channel at NLO in QCD with NNLL soft gluon corrections [63,64], and to NLO in QCD for the t-and s-channels [63,64]. Top-quark pair production with an additional W or Z boson was calculated using MADGRAPH5_aMC@NLO 2.2.2 [65] at NLO in the ME calculations, while fragmentation and hadronization were simulated with PYTHIA 8.186. The underlying-event tune A14 was used with the NNPDF2.3LO PDF set, and the cross sections were normalized using NLO predictions [66,67].
Small contributions from Higgs boson events produced by gluon-gluon fusion and vector-boson fusion were modeled using POWHEG-BOX v2 with the NNPDF3.0NNLO PDF and showered using PYTHIA 8.186. Contributions from the associated production of a Higgs boson with a vector boson and from a Higgs boson in association with two top quarks were simulated using PYTHIA 8.186 and MADGRAPH5_aMC@NLO, respectively. All Higgs boson samples were normalized to cross sections from Ref. [68].
For all samples showered with PYTHIA 8, EVTGEN 1.2.0 [69] was used to simulate the decays of bottom and charmed hadrons.

B. Simulated signal samples
Simulated signal samples were generated using MADGRAPH5_aMC@NLO 2.6.2 interfaced to PYTHIA 8.186 with the A14 tune for the PS modeling, hadronization, and underlying event. The ME calculation was performed at tree level and includes the emission of up to two additional partons. The NNPDF2.3LO PDF set was used. The ME-PS matching used the CKKW-L [70] prescription, with a matching scale set to one quarter of the mass of the pair of produced particles. Signal cross sections were calculated to next-to-leading order in the strong coupling constant, adding the resummation of soft gluon emission at next-to-leading-logarithm accuracy (NLO þ NLL) using the PDF4LHC15mc PDF set [71]. The nominal cross section and its uncertainty were taken from an envelope of cross section predictions using different PDF sets and factorization and renormalization scales, as described in Refs. [71][72][73][74][75].
The masses of all charginos and neutralinos apart from theχ 0 1 , were set to 2.5 TeV, thus leaving a single kinematically allowed decay:τ →χ 0 1 τ. The stau-left and stau-right were combined and have the same mass, which was varied between 100 and 440 GeV and no mixing was assumed between the gauge and mass eigenstates. The mass of the binolikeχ 0 1 was varied in the range 1-200 GeV. Reference points withτ masses of 120 GeV, 280 GeV, and aχ 0 1 mass of 1 GeV are used throughout this paper to illustrate typical features of the SUSY models to which this analysis is sensitive. The theoretical cross section used at NLO þ NLL was 140 (50) fb withτ LτL (τ RτR ) of 120 GeV, and 5.8 (2.2) fb withτ LτL (τ RτR ) of 280 GeV.

IV. EVENT RECONSTRUCTION
Events with at least one reconstructed primary vertex [76] are selected. A primary vertex must have at least two associated charged-particle tracks with transverse momentum p T > 500 MeV and be consistent with the collision region. In events with multiple primary vertices, the one with the largest P p 2 T of the associated tracks is chosen. Jets are reconstructed from three-dimensional calorimeter energy clusters [77] using the anti-k t algorithm [78,79] with a radius parameter of 0.4. Jet energies are corrected for detector inhomogeneities, the noncompensating response of the calorimeter, and the impact of pileup, using factors derived from test beam and pp collision data, and from a detailed GEANT4 detector simulation [80]. The impact of pileup is accounted for using a technique, based on jet areas, that provides an event-by-event and jet-by-jet correction [81]. Jets that are likely to have originated from pileup are not considered in this analysis [82]. The efficiency of this pileup-rejection selection is approximately 92%. Jets are required to have p T > 20 GeV and jηj < 2.8. Events containing jets that are likely to have arisen from detector noise or cosmic rays are removed.
Jets containing b-hadrons (b-jets) are identified using the MV2c10 algorithm [83], based on a multivariate discriminant making use of track impact parameters and reconstructed secondary vertices [84]. Candidate b-jets are required to have p T > 20 GeV and jηj < 2.5. A working point with an average b-tagging efficiency of 77% for simulated tt events is used. This working point corresponds to a c-jet and light-jet rejection of 4.9 and 110, respectively [85,86].
Electron candidates are reconstructed by matching clusters in the electromagnetic calorimeter with chargedparticle tracks in the inner detector. Electrons are required to have p T > 17 GeV and jηj < 2.47, and to satisfy the "loose" working point according to a likelihood-based identification [87]. Muon candidates are reconstructed from MS tracks matching ID tracks. Muons are required to have p T > 14 GeV and jηj < 2.7 and fulfill the "medium" quality criteria of Ref. [88]. Events containing a muon candidate with a poorly measured charge-to-momentum ratio (σðq=pÞ=jq=pj > 0.4) are rejected.
Electrons and muons are required to satisfy isolation criteria to reduce the number of jets misidentified as charged leptons. The scalar sum of the p T of tracks inside a variable-size cone around the lepton (excluding its own track), must be less than 15% of the lepton p T . The track isolation cone size 2 for electrons (muons) ΔR is given by the minimum of ΔR ¼ 10 GeV=p T and ΔR ¼ 0.2 (0.3). In addition, for electrons (muons) the sum of the transverse energy of the calorimeter energy clusters in a cone of ΔR ¼ 0.2 around the lepton, excluding the energy from the lepton itself, must be less than 20% (30%) of the lepton p T . For electrons with p T > 200 GeV these isolation requirements are not applied, and instead an upper limit of maxð0.015 × p T ½GeV; 3.5 GeVÞ is placed on the transverse energy of the calorimeter energy clusters in a cone of ΔR ¼ 0.2 around the electron.
The efficiencies for electrons and muons to satisfy the reconstruction, identification, and isolation criteria are measured in samples of leptonic Z and J=ψ decays, and corrections are applied to the simulated samples to reproduce the efficiencies in data.
The reconstruction of hadronically decaying τ-leptons is based on information from tracking in the ID and threedimensional clusters of energy in the electromagnetic and hadronic calorimeters. The τ-lepton reconstruction algorithm is seeded by jets reconstructed as described above but with p T > 10 GeV and jηj < 2.5. The reconstructed energies of the hadronically decaying τ-lepton candidates are corrected to the τ-lepton energy scale, which is calibrated using simulation and in situ measurements using Z → ττ decays [89]. The τ-neutrino from the hadronic τ-lepton decay is not taken into account in the reconstruction and calibration of the τ-lepton energy and momentum. Hadronic τ-lepton decay candidates are required to have one or three associated charged-particle tracks (prongs) and the total electric charge of those tracks must be AE1 times the electron charge. To improve the discrimination between hadronically decaying τ-leptons and jets a multivariate algorithm is used [89]. The τ-lepton identification algorithm is based on a boosted decision tree (BDT) method. The BDT algorithm uses various track and cluster variables as input to discriminate between τ-leptons and jets. For 1-prong (3-prong) τ-lepton candidates, the signal efficiencies are 75% (60%) and 60% (45%) for the "medium" and "tight" working points, respectively. In the following, τ-lepton candidates are required to satisfy the "medium" identification criteria for jet discrimination (medium τ-lepton candidates), unless otherwise stated. For electron discrimination, an additional BDT is trained to discriminate τ-leptons from electrons and a selection is applied on the output of the discriminant. This requirement has about 95% efficiency for selecting τ-leptons, and a background electron rejection factor from 10 to 50 depending on the η. The τ-lepton candidates are required to have p T > 20 GeV and jηj < 2.5, excluding the transition region between the barrel and end cap calorimeters (1.37 < jηj < 1.52).
The simulation is corrected for differences between the efficiencies of τ identification in data and simulation at both trigger and reconstruction level. For hadronically decaying τ-leptons originating from prompt gauge-boson decays, the corrections are calculated with a tag-and-probe method in a sample of Z → ττ events where one τ-lepton decays hadronically and the other leptonically into a muon and two neutrinos [90]. Small differences between data and simulation for electron, muon, jet identification, and trigger efficiencies are corrected with scale factors derived in control regions. For the trigger corrections, dedicated regions enriched in tt, Z, and W events are used.
The measurement of the missing transverse momentum vector, p miss T , and its magnitude, E miss T , is based on the negative vectorial sum of the p T of all identified jets, τ-lepton candidates, electrons, muons, photons [87], and an additional soft term. The soft term is constructed from all high-quality tracks that are associated with the primary vertex but not with any identified particle or jet. In this way, the missing transverse momentum is adjusted for the best calibration of the jets and the other identified particles, while maintaining pileup independence in the soft term [91].
Possible double counting of reconstructed objects is resolved in the following order. The τ-lepton candidates close to electron or muon candidates (ΔR < 0.2) are removed, as are electrons that share a track with a muon. For electrons close to a jet (ΔR < 0.4), the electron is removed, except when ΔR < 0.2 and the jet is not b-tagged, in which case the jet is removed. For a muon close to a jet (ΔR < 0.4), the muon is removed unless the jet has less than three tracks associated with it and is within ΔR < 0.2. In the latter case, the jet is removed. Any remaining jet within ΔR ¼ 0.2 of a τ-lepton candidate is removed.

V. EVENT SELECTION
The events are required to have exactly two medium τ-lepton candidates with opposite-sign electric charge (OS) and to have passed either an asymmetric di-τ trigger or a combined di-τ þ E miss T trigger [92]. In events selected by the di-τ þ E miss T trigger, the offline reconstructed E miss T must be larger than 150 GeV, to achieve full trigger efficiency. In 2015-2017 data the efficiency of the asymmetric di-τ (di-τ þ E miss T ) trigger is 75-80% for events with the leading τ-lepton candidate p T > 95 (50) GeV and the subleading , where y is the rapidity and ϕ the azimuthal angle.
τ-lepton candidate p T > 60ð40Þ GeV. The efficiency is evaluated for events with correctly identified τ-leptons and the p T thresholds are quoted at reconstruction level. As the luminosity increased in later years, the trigger selection was tightened; in 2018 the same efficiency was reached for the leading τ-lepton candidate p T > 95ð75Þ GeV and the subleading τ-lepton candidate p T > 75ð40Þ GeV for the asymmetric di-τ (di-τ þ E miss T ) trigger. Events with an additional third medium τ-lepton or light lepton are rejected.
The reconstructed invariant mass of the visible decay products of the two leading τ-lepton candidates, mðτ 1 ; τ 2 Þ, must be larger than 120 GeV to remove τ-leptons originating from decays of low-mass resonances and to suppress contributions from Z þ jets and Higgs boson events (Z=H-veto).
In order to further discriminate between SUSY signal events and SM background processes, additional requirements are applied to define the signal region (SR) selections. To reject events from SM processes containing a top quark, selected events must not contain any b-tagged jets (b-jet veto). A lower bound on the stransverse mass m T2 [93,94] is imposed to reduce contributions from tt and WW events. The m T2 variable is defined as: where p T;τ 1 and p T;τ 2 are the transverse momenta of the two τ-lepton candidates, and the transverse momentum vector of one of the invisible particles, q T , is chosen to minimize the larger of the two transverse masses m T;τ 1 and m T;τ 2 . The latter masses are defined by For tt and WW events, in which two W bosons decay leptonically and p miss T is the sum of the transverse momenta of the two neutrinos, the m T2 distribution has a kinematic endpoint at the W boson mass. For large mass differences between the staus and the lightest neutralino, the m T2 distribution for signal events extends significantly beyond this endpoint.
The SRs were optimized for stau discovery by varying the kinematic selection criteria. Events are selected by the asymmetric di-τ trigger to cover the low stau mass region (SR-lowMass) and by the di-τ þ E miss T trigger to cover the high stau mass region (SR-highMass), as shown in Table I. Values of E miss T < 150 GeV (E miss T > 150 GeV) are required for SR-lowMass (SR-highMass) to keep the two SRs orthogonal. Both (at least one) of the τ-lepton candidates must satisfy the tight identification criteria for jet discrimination ('tight' τ-lepton candidate) for SR-lowMass (SR-highMass). Values of E miss T > 75 GeV are required for SR-lowMass to increase signal sensitivity.

VI. STANDARD MODEL BACKGROUND ESTIMATION
The main SM processes contributing to the selected final states are multijet, W þ jets and diboson production. Background events may contain a combination of "real" τ-leptons, defined as correctly identified prompt τ-leptons, or "fake" τ-leptons, which can originate from a misidentified quark or gluon jet, an electron, or a muon.
In multijet events in the SRs, nearly all τ-lepton candidates are misidentified jets. The multijet contribution in the SRs is estimated from data, as described in Sec. VI A. The contribution arising from heavy-flavor multijet events containing a real τ-lepton from the heavy-flavor quark decay is included in the multijet estimation. The contribution of W þ jets events, which contain one real τ-lepton from the W boson decay and one or more misidentified jets, is estimated from MC simulation and normalized to data in a dedicated control region (CR), as described in Sec. VI B.
Diboson production contributes mainly events containing real τ-leptons originating from WW and ZZ decaying into a ττνν final state. Additional SM backgrounds arise from Z þ jets production, or events that contain a top quark or a top-quark pair in association with jets or additional W or Z bosons (collectively referred to as top-quark background in the following). The contribution from real τ-leptons exceeds 90% in Z þ jets and diboson production, and ranges from 45% to 75% in backgrounds containing top quarks according to the MC simulation. The contribution of fake τ-leptons from heavy-flavor decays in jets is found to be negligible in MC simulation. To estimate the irreducible background, which includes diboson, Z þ jets and top-quark events, only MC simulated samples are used and validated in dedicated validation regions (VRs), as described in Sec. VI C. A simultaneous fit of the event yields in the control regions, based on the profile likelihood method, is performed to normalize the multijet and W þ jets background estimates and propagate systematic uncertainties, as described in Sec. VIII. The sources of systematic uncertainty in the background estimates are described in Sec. VII. Results showing the level of agreement in the validation regions after the fit are labeled "postfit" and results shown in the control regions prior to the fit are labeled "prefit."

A. Multijet background estimation
One of the dominant backgrounds in the SRs originates from jets misidentified as τ-leptons in multijet production. It accounts for 44% (30%) of the total SM contribution in SR-lowMass (SR-highMass). This contribution is estimated from data using the so-called ABCD method. All regions used for the ABCD method are schematically drawn in Fig. 2. Four exclusive regions, labeled as A, B, C, and D, are defined in a two-dimensional plane as a function of two (or more) discriminating variables that are largely uncorrelated. The ratio of events in the regions C and B is then equal to that in the regions D and A. The number of multijet events in region D, N D , can therefore be calculated from the multijet events in region A, N A , multiplied by the transfer factor T ¼ N C =N B , where N C (N B ) is the number of multijet events in region C (B). The region D corresponds to the SR defined in Sec. V, whereas the regions A, B, and C are control regions defined accordingly. In the following, the regions A, B, C, D are labeled as CR-A, CR-B, CR-C and SR-lowMass (or SR-highMass), respectively. The ABCD method only provides a first-order estimate of multijet background, the normalization and uncertainty being then modified by a combined fit to CR-A (lowMass) and CR-A (highMass) described in Sec. VIII, in which the transfer factor is used as an input and fixed. The above two-fold method has the advantages of taking the possibility of contamination from signal in multijet CR-A into account in the signal exclusion fit described in Sec. X, as well as considering the correlation of systematic uncertainties among control regions and among background processes.
The definition of the regions used in the ABCD method for the multijet estimation is given in Table II; only those requirements that differ between the CRs/VRs and the SRs are listed. In all of the regions, the di-τ þ E miss T or asymmetric di-τ triggers described in Sec. V are used. CR-A and CR-B include the events with two "loose" τ-leptons, either same sign (SS) or OS, and veto the events with two medium OS τ-leptons to remain orthogonal to the SR and reduce the potential contamination by signal events. CR-C and CR-B have events with lower m T2 and E miss T ; no ΔRðτ 1 ; τ 2 Þ requirement is applied. Furthermore, two sets of VRs, VR-E and VR-F, are defined corresponding to each SR. The validation regions are used to validate the extrapolation of the ABCD estimation to the SRs and to estimate the systematic uncertainty from the residual correlation between the τ-lepton identification and charge requirements and the kinematic variable m T2 .
The number of multijet events in each control region and validation region is estimated from data after subtraction of other SM contributions estimated from MC simulation. In CR-B and VR-E, respectively, around 96% and 90% (75% and 79%) of the events are from multijet production in the lowMass (highMass) region. For CR-A and CR-C, respectively, the multijet purity is 74% and 57% (58% and 53%) FIG. 2. Illustration of the ABCD method for the multijet background determination for SR-lowMass (left) and SR-highMass (right). The control regions A, B, C, and signal region D for the ABCD method described in the text (labeled as CR-A, CR-B, CR-C and SR-lowMass/SR-highMass) are drawn as white boxes. Shown shaded and labeled VR are the regions E and F, which are used to validate the ABCD method and to estimate the systematic uncertainty. The transfer factor T used in the ABCD method is the ratio of number of multijet events in the regions C and B.
Good agreement between data and the estimated SM background is found for the m T2 and E miss T distributions in the validation regions, as shown in Fig. 3.
The contamination from signal in multijet CR-A is defined as the ratio of the number of signal events to the sum of the number of signal events and SM background processes. It ranges from 0.4% (1.2%) to 9.4% (21.4%) for SR-lowMass (SR-highMass) and is taken into account in the fit.
The multijet background estimation is also validated using a different method, the fake-factor method. Fake factors (FF) are derived in a region enriched in multijet events consisting of events with two τ-lepton candidates with the same electric charge. Both τ-lepton candidates are used to compute the fake factor as the ratio to the number of leading or subleading τ-lepton candidates satisfying all of the nominal signal identification criteria over the number of τ-lepton candidates failing the requirement on the jetveto BDT, and are parametrized by the p T , η and N prong of the τ-lepton candidates. An observed dependence on the BDT output score of the second τ-lepton candidate in the event is also taken into account. The expected number of multijet events entering a selection is computed by applying the fake factors to the yields in a set of sideband regions where either the leading, subleading or both τ-lepton candidates fail the identification requirements. The predicted multijet event yields from the ABCD method and the FF method in both SRs and the multijet VRs agree within the statistical and systematic uncertainties.

B. W + jets background estimation
The production of W þ jets events with at least one misidentified τ-lepton accounts for about 25% of the expected SM background in the two SRs. A dedicated control region (WCR) is used to normalize the W þ jets MC estimate to data and another region is used to validate the W þ jets estimate (W validation region, WVR). To suppress contamination by multijet events, the WCR is enriched in events where the W boson decays leptonically into a muon and a neutrino. Events are selected with a single-muon trigger using the lowest unprescaled p T thresholds available. Events containing exactly one muon and one candidate τ-lepton with opposite electric charge are selected. The muon is required to have p T > 50 GeV. The τ-lepton candidate must satisfy the medium τ-lepton identification criteria and is required to have p T > 60 GeV.
The contribution from events with top quarks is suppressed by rejecting events containing b-tagged jets, and rejecting events which are kinematically compatible with tt production (top-tagged) through the use of the contransverse mass variable m CT [95]. The definitions of the WCR and WVR are given in Table III. The transverse mass of the μ þ E miss T system, m T;μ is used to reduce the contribution from Z þ jets, top-quark, and diboson events. To further suppress multijet and Z þ jets events, E miss T and ΔRðτ; μÞ cuts are applied. The invariant mass of the muon and τ-lepton, mðμ; τÞ, and m T;τ þ m T;μ are used to improve the W þ jets purity. Events in the WCR (WVR) are selected by requiring low (high) m T2 .
The multijet contribution in the WCR (WVR) is estimated using the so-called OS-SS method by counting the number of events in data satisfying the same requirements as for the WCR (WVR) but with the electric charge of the two leptons having the same sign. Events from SM TABLE II. Definition of the regions used in the ABCD method for the multijet estimation in the SRs. Only those requirements that differ between the CRs/VRs and the SRs are listed.

CR-B (highMass) CR-C (highMass)
≥2 loose τ 2 medium τ (OS) <2 medium τ (OS) ≥1 tight τ no ΔRðτ 1 ; τ 2 Þ cut no ΔRðτ 1 ; τ 2 Þ cut 50 < E miss T < 100 GeV 50 < E miss T < 100 GeV 30 < m T2 < 50 GeV 30 < m T2 < 50 GeV processes other than multijet production are subtracted from the data counts in the SS region using MC simulation. The OS-SS method relies on the fact that in the multijet background the ratio of SS to OS events is close to unity, while it is around one-seventh for W þ jets production. The latter is dominated by gu=gd-initiated processes that often give rise to a jet originating from a quark, the charge of which is anticorrelated with the W boson charge. Based on studies with simulated samples, a systematic uncertainty of 100% is assigned to the multijet estimate in the WCR. The purity of the W þ jets selection is around 79% (69%) in the WCR (WVR). The prefit m T2 distribution in the WCR is shown in Fig. 4, and good agreement, both for the normalization and shape, between data and SM predictions is observed. The contamination from signal in the WCR and WVR is negligible.

C. Irreducible background estimation
Irreducible SM backgrounds arise mainly from tt, single top quark, tt þ V, Z þ jets, and multiboson (diboson (WW, WZ and ZZ), triboson (VVV)) and Higgs boson. They are estimated with MC simulation. Other SM backgrounds are found to be negligible.
The total contribution from tt, single top quarks, tt þ V and Z þ jets amounts to about 8% and 20% of the total background in SR-lowMass and SR-highMass, respectively. The diboson background accounts for 23-25% of the

WCR WVR
1 medium τ and 1 isolated μ (OS) single-muon trigger p T ðτÞ > 60 GeV, p T ðμÞ > 50 GeV E miss T > 60 GeV b-jet veto and top-tagged events veto mðμ; τÞ > 70 GeV 1 < ΔRðμ; τÞ < 3.5 50 < m T;μ < 150 GeV m T;μ þ m T;τ > 250 GeV 30 < m T2 < 70 GeV m T2 > 70 GeV total SM contribution in the SRs and mainly arises from WW → τντν and ZZ → ττνν events, in which more than 96% of the contribution is from events with two real τ-leptons according to the MC simulation. The MC estimates are validated in regions enriched in Z þ jets, top-quark, and multiboson events. For these validation regions, events are required to pass either the combined di-τ þ E miss T trigger or the asymmetric di-τ trigger with the same offline threshold as described in Table IV. Events are required to have at least two τ-lepton candidates satisfying the medium τ-lepton identification criteria with opposite electric charge, and at least one τ-lepton candidate must satisfy the "tight" τ-lepton identification criteria in order to be close to the SRs. In the top-quark validation regions (TVR), to increase the contribution from top-quark events, events must satisfy the requirement ΔRðτ 1 ; τ 2 Þ > 1.2 and contain at least one b-tagged jet with p T > 20 GeV. In order to be close to the SRs, m T2 > 60 GeV is required. In the Z þ jets and multiboson validation regions (ZVR, VVVR), in order to suppress top-quark backgrounds, events containing b-tagged jets are vetoed. To further enhance the purity of Z þ jets events, ΔRðτ 1 ; τ 2 Þ, mðτ 1 ; τ 2 Þ and m T2 requirements are applied. Further, mðτ 1 ; τ 2 Þ, m T;τ1 þ m T;τ2 and m T2 cuts are applied in order to enhance multiboson purity. The ZVRs, TVRs and VVVRs requirements are summarized in Table IV. The data event yields and the SM predictions in the W þ jets (WVR), top quark (TVR), Z þ jets (ZVR) and multiboson (VVVR) validation regions are shown in Fig. 5. The data and SM prediction in each validation region agree within the uncertainties. The purity of the selection in Z þ jets and tt (multiboson) events is 83%-96% (47%-71%) in the respective validation regions.

VII. SYSTEMATIC UNCERTAINTIES
Systematic uncertainties have an impact on the estimates of the background and signal event yields in the control and signal regions. Uncertainties arising from experimental and theoretical sources are estimated.
The main sources of experimental systematic uncertainty in the SM background estimates include τ-lepton and jet energy scale and resolution, τ-lepton identification, pileup, and uncertainties related to the modeling of E miss T in the simulation. The uncertainties in the energy and momentum scale of each of the objects entering the E miss T calculation are estimated, as well as the uncertainties in the soft-term resolution and scale. A pileup reweighting procedure is applied to simulation to match the distribution of the number of reconstructed vertices observed in data [96]. The corresponding uncertainty is derived by a reweighting in which <μ> is varied by AE4%. The main contributions to experimental systematic uncertainties in the SRs are from the τ-lepton identification and energy scale [89], and jet energy scale and resolution. Systematic uncertainties associated with the τ-lepton triggers are also included GeV τp T and E miss T cuts described in Sec. V with the τ-lepton identification. Other contributions are less than 3%.
Theoretical uncertainties affecting the SHERPA MC event generator predictions, for the W þ jets, Z þ jets and multiboson samples, are estimated by varying the PDF sets following the PDF4LHC recommendations [97] as well as the seven variations of the QCD renormalization and factorization scales in the matrix element by factors of 2 and 0.5 avoiding variations in opposite direction [98]. For tt production, uncertainties in the parton shower simulation, modeling of initial-and final-state radiation, and associated with the choice of event generator are considered. The theory uncertainties are implemented as a separate nuisance parameter for each scale variation and process; they are not constrained. The theory uncertainty in W þ jets (multiboson) production is mainly due to the QCD renormalization scale variation, which amounts to 2%-3% (5%-6%) compared to the total background yield.
The following sources of uncertainty are considered for the ABCD method used to determine the multijet background: the correlation among the τ-lepton identification, the charge requirement, and the kinematic variable m T2 , the limited number of events in the CRs, and the subtraction of other SM backgrounds. The systematic uncertainty in the correlation is estimated by comparing the transfer factor from CR-B to CR-C to that of VR-E to VR-F. The combined experimental systematic uncertainty and theory uncertainty in the nonmultijet background subtraction in the control regions is estimated by considering the systematic uncertainty of the MC estimates of the nonmultijet background in these regions. These uncertainties are negligible due to high multijet purity in the CRs. The statistical uncertainty of the event yields in the control regions is propagated to the signal regions as a systematic uncertainty.
The systematic uncertainties of the background estimates in the SRs are summarized in Table V. The dominant uncertainties in the SRs are the statistical uncertainty of the MC predictions (11%-21%), τ-lepton identification and energy scale (10%-19%), and multijet background normalization (8%-12%).
The total uncertainty in the signal yields for the SUSY reference points defined in Sec. III is about 17%-31%. The dominant uncertainties in the SRs are from τ-lepton The SM backgrounds other than multijet production are estimated from MC simulation. The multijet contribution is negligible and is estimated from data using the ABCD method, using CRs obtained with the same technique used for the SRs, and described in Sec. VI A. The hatched bands represent the combined statistical and systematic uncertainties of the total SM background. The background-only fit to data is used, described in Sec. VIII. For illustration, the distributions from the SUSY reference points (defined in Sec. III) are also shown as dashed lines. The lower panel shows the ratio of data to the SM background estimate. identification and the τ-lepton energy scale (14%-29%) and the statistical uncertainty of the signal MC predictions (6%-10%). The cross section uncertainty is taken into account as the main source of the signal modeling theoretical uncertainty; and it varies from 2% to 3% for the SUSY models considered.

VIII. STATISTICAL ANALYSIS
The statistical interpretation of the results is performed using the profile likelihood method implemented in the HistFitter framework [99]. The likelihood function is a product of the probability density functions for every region contributing to the fit. The number of events in a given CR or SR is described using a Poisson distribution, the mean of which is the sum of the expected contributions from all background and signal sources. The systematic uncertainties in the expected event yields are included as nuisance parameters and are assumed to follow a Gaussian distribution with a width determined by the size of the uncertainty. Correlations between control and signal regions, and among background processes are taken into account through common nuisance parameters. The fit parameters are determined by maximizing the product of the Poisson probability functions and the Gaussian probability constraints for the nuisance parameters.
Three types of fits are performed for the combined SR-lowMass and SR-highMass regions.
(i) The background-only fit uses as input the number of observed events and expected SM contributions in the multijet CR-A and WCR as well as the transfer factors, which relate the number of multijet or W þ jets events in the control regions to the number in the signal regions. The free parameters in the fit are the normalizations of the W þ jets and TABLE VI. Observed and expected numbers of events in the control and signal regions. The expected event yields of SM processes are given after the background-only fit described in Sec. VIII. The entries marked as "-" are negligible. The uncertainties correspond to the sum in quadrature of statistical and systematic uncertainties. The correlation of systematic uncertainties among control regions and among background processes is fully taken into account.  (ii) A model-independent limit fit combines the data event yield in a given SR with the SM background estimate and its uncertainties to test whether any non-SM signal contributes to the SR. The significance of a possible excess of observed events over the SM prediction is quantified by the one-sided probability, pðsignal ¼ 0Þ denoted by p 0 , of the background alone to fluctuate to the observed number of events or higher using the asymptotic formula described in Ref. [100]. The background yields and uncertainties are taken from the background-only fit results. The presence of a non-SM signal would manifest itself as a small p 0 value. (iii) In the model-dependent limit fit the SUSY signal is allowed to populate both the signal regions and the control regions, and it is scaled by a freely floating signal normalization factor. The background normalization factors are also determined simultaneously in the fit. A SUSY model with a specific set of sparticle masses is rejected if 95% C.L. upper limit on the signal normalization factor obtained in this fit is smaller than unity.

IX. RESULTS
The observed number of events in each control, and signal region and the expected contributions from SM processes are given in Table VI. The contributions of multijet and W þ jets events are scaled with the normalization factors obtained from the background-only fit described in Sec. VIII. The multijet normalization in the SR is 1.03 times the prediction from the ABCD method and has an uncertainty of around 28%, due to the small number of observed events in the multijet CR-A. The W þ jets normalization is 0.91 AE 0.12. The m T2 distributions are shown in Fig. 6 for data, expected SM backgrounds, and the SUSY reference points defined in Sec. III. In both signal regions, observations and background predictions are found to be compatible within uncertainties.
Individual model-independent fits of the SR-lowMass and SR-highMass are used to derive the one-sided p 0values. In addition, the observed and expected 95% C.L. upper limits on the visible non-SM cross section (σ 95 vis ) are derived using the same model-independent fits. The σ 95 vis is defined as the product of acceptance, reconstruction efficiency and production cross section, obtained normalizing the limit on the non-SM yield in the SR by the integrated luminosity of the data sample. These results are shown in Table VII. All limits are calculated using the CL s prescription [101].

X. INTERPRETATION
In the absence of a significant excess over the expected SM background, the observed and expected numbers of events in the signal regions are used to set exclusion limits at 95% C.L. using the model-dependent limit fit. The exclusion limits for the combined SR-lowMass and SR-highMass for the simplified models described in Sec. III are shown in Fig. 7. The solid (dashed) lines show the observed (expected) exclusion contours. The band  around the expected limit shows the AE1σ variations, including all uncertainties except theoretical uncertainties in the signal cross section. The dotted lines around the observed limit indicate the sensitivity to AE1σ variations of the theoretical uncertainties in the signal cross section. Stau masses from 120 GeV to 390 GeV are excluded for a massless lightest neutralino in the scenario of combined stau-left and stau-rightτ þ R;Lτ − R;L production. For stau-left pair production only,τ þ Lτ − L , the exclusion region extends from 155 GeV to 310 GeV. While the stau-left pairs have a higher production cross section, the stau-right pairs have a higher efficiency times acceptance due to kinematic differences in the resulting decay products.

ATLAS
These limits extend significantly beyond previous results [25,27,102,103] in the highτ mass region.
Uncertainties in each background from scale variations are fully correlated across regions and bins, and uncorrelated between processes. In some cases this may result in uncertainties canceling out, while the higher-order corrections may not cancel out. A different fit with scale variations uncorrelated in all bins results in a less than 1% change in the limits on the excluded cross section near the edge of exclusion for combined left and rightττ production.

XI. CONCLUSION
Searches for stau-pair production in events with at least two hadronically decaying τ-leptons and missing transverse momentum were performed using 139 fb −1 of pp collision data at ffiffi ffi s p ¼ 13 TeV recorded with the ATLAS detector at the LHC. Agreement between data and SM predictions is observed in two optimized signal regions. The results are used to set limits on the visible cross section for events beyond the Standard Model in each signal region. Exclusion limits are placed on parameters of simplified electroweak supersymmetry models in scenarios of staupair production. Stau masses from 120 GeV to 390 GeV are excluded for a massless lightest neutralino in the scenario of direct production of stau pairs, with each stau decaying into the lightest neutralino and one τ-lepton. These limits extend significantly beyond previous results by the ATLAS and CMS experiments in the highτ mass region. lepsusy/Welcome.html.