Search for long-lived particles decaying to leptons with large impact parameter in proton-proton collisions at $\sqrt{s}$ = 13 TeV

A search for new long-lived particles decaying to leptons using proton-proton collision data produced by the CERN LHC at $\sqrt{s}$ = 13 TeV is presented. Events are selected with two leptons (an electron and a muon, two electrons, or two muons) that both have transverse impact parameter values between 0.01 and 10 cm and are not required to form a common vertex. Data used for the analysis were collected with the CMS detector in 2016, 2017, and 2018, and correspond to an integrated luminosity of 118 (113) fb$^{-1}$ in the ee channel (e$\mu$ and $\mu\mu$ channels). The search is designed to be sensitive to a wide range of models with displaced e$\mu$, ee, and $\mu\mu$ final states. The results constrain several well-motivated models involving new long-lived particles that decay to displaced leptons. For some areas of the available phase space, these are the most stringent constraints to date.


Introduction
To date, no evidence of particles beyond the standard model (BSM) has been found by any experiment, including those at the CERN LHC. However, the vast majority of LHC searches assume the lifetimes of new particles are short enough that their decay products are prompt, i.e., their decay products are consistent with originating at the primary proton-proton (pp) interaction. Search strategies are rarely optimized for particles with measurable lifetimes whose decays produce detector signatures that are displaced from the primary pp interaction. Therefore, new phenomena with displaced signatures can escape such searches. Particles in BSM scenarios can have long lifetimes for the same reasons as standard model (SM) particles, namely, small couplings between the long-lived particles (LLPs) and lighter states, approximate symmetries, heavy intermediate states, or limited phase space availability for decays.
The CMS Collaboration has previously performed a search for signatures with one displaced electron and one displaced muon, using pp collision data recorded at √ s = 8 TeV and corresponding to an integrated luminosity of 19.7 fb −1 [26]. Both this previous analysis and the one described in this paper, which uses data recorded at √ s = 13 TeV, are optimized to the phase space just beyond the sensitivity of prompt searches but with smaller displacements than other searches for long-lived BSM signatures. As a result, this search is sensitive to LLPs with proper decay lengths (cτ 0 ) between approximately 10 -3 and 10 3 cm, where c is the speed of light and τ 0 is the proper lifetime. These two analyses are unique in that they allow but do not require the displaced final-state particles to originate from a common vertex. Such a vertex is often required in other searches under the assumption that an LLP will decay to multiple leptons. In contrast, here we perform an inclusive search that is sensitive to one LLP that decays to multiple leptons, two LLPs that each decay to at least one lepton, and any other topology whose final state includes at least two displaced leptons. The search described in this paper uses data taken with the CMS detector during 2016, 2017, and 2018. We conduct a search for events in which there are one electron and one muon, two electrons, or two muons in the final state, where both of the leptons are displaced from the beam axis. With respect to the previous search, this search is based on about a factor of 6 larger integrated luminosity, is performed at a higher √ s, and adds two same-flavor channels corresponding to the ee and µµ final states.
This search is designed to be model independent and to be sensitive to as many event topologies as possible. Consequently, the event selection focuses exclusively on a displaced, isolated dilepton signature and does not try to identify signal events using hadronic activity or missing transverse momentum from undetected particles. In this way, we retain sensitivity to any model that can produce leptons with displacements on the order of 0.01 to 10 cm and with sufficiently high momenta, regardless of whether these leptons are accompanied by jets, missing transverse momentum, or other kinematic features.
We interpret the search results in the context of several models that produce final states with displaced leptons, the Feynman diagrams for which are shown in Fig. 1. The first model introduces R-parity-violating (RPV) terms in the superpotential of the minimally supersymmet-ric (SUSY) SM [27], allowing the lightest SUSY particle (LSP) to decay to SM particles. Only lepton-number-violating operators are considered because of constraints from measurements of the lifetime of the proton [28]. With sufficiently small couplings for these operators, the LSP has a long enough lifetime that its decay products are measurably displaced from the pp interaction region. We focus on the case in which the LSP is the top squark ( t), the superpartner of the top quark. At the LHC, the top squark would be dominantly pair produced, and we consider its decay through an RPV vertex to a d (b) quark and a charged lepton via t → d ( t → b ). With respect to the previous analysis [26], we have added the t → d decay to facilitate reinterpreting the results in a wider range of BSM scenarios, although we find that this decay mode produces similar results to those from the t → b decay. We expect that the t → s decay mode will also produce similar results, although we did not explore it. For simplicity, we assume lepton universality in the top squark decay vertex, so that the branching fraction to any lepton flavor, namely, an electron, muon, or tau lepton, is equal to one-third. We also interpret the results with a gauge-mediated SUSY breaking (GMSB) model in which the next-to-lightest SUSY particle (NLSP) is long lived because of its small gravitational coupling to the LSP gravitino G, which is nearly massless [29]. In this model, the NLSP is the superpartner of a lepton (slepton ). We consider selectrons e, smuons µ , and staus τ separately, as well as together, as co-NLSPs. The sleptons would be pair produced at the LHC and each would decay to a lepton (e, µ, or τ) of the same flavor and a gravitino. In addition, we consider a model that produces BSM Higgs bosons (H) with a mass of 125 GeV through gluon-gluon fusion [30]. The H decays to two long-lived scalars S, each of which decays to two oppositely charged and same-flavor leptons, and the probabilities of the lepton pair being ee or µµ are assumed to be the same. For the scenarios where the long-lived top squarks and sleptons decay to tau leptons, events in which the tau leptons both subsequently decay into electrons or muons are considered.
Tabulated results are provided in the HEPData record for this analysis [31].

The CMS detector
The central feature of the CMS apparatus is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter, each composed of a barrel and two endcap sections. The electron momentum is estimated by combining the energy measurement in the ECAL with the momentum measurement in the tracker. Forward calorimeters extend the pseudorapidity η coverage provided by the barrel and endcap detectors. Muons are measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid using three technologies: drift tubes (DTs) in the barrel, cathode strip chambers (CSCs) in the endcaps, and resistive-plate chambers in both the barrel and the endcaps. Each of the four muon detector planes provides reconstructed hits on several detection layers, which are combined into local track segments, forming the basis of muon reconstruction inside the muon system.
The silicon tracker measures charged particles with |η| < 3.0. During the 2016 LHC run, the silicon tracker consisted of 1440 silicon pixel and 15 148 silicon strip detector modules. The pixel detector was then upgraded, such that in the 2017 and 2018 LHC runs, it consisted of 1856 pixel modules. After the upgrade, the number of pixel layers increased from three, with radii between 4.4 and 10.2 cm from the interaction region, to four, with radii between 2.9 and 16.0 cm [32]. In the 2017 and 2018 LHC runs, for nonisolated particles within the transverse momentum range 1 < p T < 10 GeV and |η| < 3.0, the average track p T resolution is 1.5%. The transverse impact parameter (d 0 ) is defined as the distance of closest approach in the transverse plane of the helical trajectory of the track with respect to the beam axis [33,34]. The sign of d 0 is given by the sign of the scalar product between the lepton momentum and the vector from the beam axis to the lepton track reference point. The track d 0 resolution improves from 75 to 20 µm as the p T increases, for nonisolated particles with 1 < p T < 10 GeV and |η| < 3.0 in 2017 and 2018. With the upgraded silicon pixel tracker, the d 0 resolution is approximately 25% better than in earlier data sets. The efficiency to reconstruct tracks as a function of |d 0 | is given in Ref. [34].
Events of interest are selected using a two-tiered trigger system. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a fixed latency of about 4 µs [35]. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage [36,37].
A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in Ref. [38].

Data and Monte Carlo simulation
The data correspond to an integrated luminosity of 118 (113) fb −1 in the ee channel (eµ and µµ channels), 16 fb −1 of which were collected in 2016 before the pixel detector upgrade. The difference in integrated luminosity in the different channels is due to the availability of the triggers, which will be described in Section 5. In addition, we use events containing muons from cosmic ray showers that were recorded with dedicated triggers for a control sample to evaluate the tracking efficiency of displaced particles, as will be described later.
In the Monte Carlo (MC) simulation of background and signal processes, minimum-bias interactions are superimposed on each event to simulate the effect of additional interactions within the same or neighboring bunch crossing (pileup). The frequency distribution of the additional interactions is adjusted to match that observed in data. The average number of pileup interactions was 23 (32) in 2016 (2017 and 2018). For the 2016 samples, the NNPDF3.0 [39] parton distribution function (PDF) set at next-to-leading order (NLO) is used, while for the samples describing the 2017 and 2018 data, the NNPDF3.1 PDF set computed at next-to-NLO order [40] is used. The PYTHIA generator is used to simulate the parton showering and hadronization for all processes. The modeling of the underlying event uses PYTHIA 8.226 [41] with the CUETP8M1 [42] tune and PYTHIA 8.230 with the CP5 tune [43] for simulated samples corresponding to the 2016 and 2017-2018 data sets, respectively. The MC-generated events are then processed through a detailed simulation of the CMS detector based on GEANT4 [44] and are reconstructed with the same algorithms used for data.
While the background is estimated using data control samples, simulated background samples are produced to perform basic checks such as comparisons of data and simulation in control regions (CRs). Samples of simulated Z+jets, W+jets, and tt production are generated at leading order (LO) using MADGRAPH5 aMC@NLO v2.2.2 (v2.4.2) for the 2016 (2017 and 2018) samples [45] and the MLM merging scheme between jets from matrix element calculations and parton showers [46]. Samples of diboson and single top quark events are simulated at NLO with POWHEG v2.0 [47][48][49][50][51]. Quantum chromodynamics (QCD) multijet events, which give rise to a background from SM events composed uniquely of jets produced through the strong interaction, are simulated with PYTHIA, selecting events that are enriched in muons. Samples of the signal process pp → t t, with the top squarks decaying via t → d or t → b , are generated using PYTHIA at LO. The top squarks can form strongly produced hadronic states called R-hadrons, which are generated with PYTHIA. In the samples used in this analysis, the interactions of the R-hadrons with matter are not simulated in GEANT4. However, such interactions are not expected to have a significant impact because these particles do not encounter a significant number of interaction lengths before decaying. Nevertheless, we study the impact of the R-hadron interactions using the "cloud model," which assumes that the top squark is surrounded by a cloud of colored, light constituents that interact during scattering [52,53], and find the effect on the signal efficiency to be negligible. The GMSB sleptons are generated at LO using MADGRAPH5 aMC@NLO v2.6.5 and PYTHIA, and the slepton decay via → G is simulated using GEANT4, which ignores information about the chirality of the slepton. Thus, we do not present results for left-and right-handed sleptons separately. The signal process pp → H → SS, S → + − is generated using POWHEG and PYTHIA at NLO.

Analysis strategy
We perform an inclusive search for displaced leptons by selecting events with wellreconstructed electrons and muons, and by rejecting background events from SM processes that produce displaced leptons, as will be described in Section 5. We use the lepton |d 0 |, which is the main discriminating variable in the analysis, to define the signal regions (SRs). Using data in a prompt CR, we correct the |d 0 | distributions in the background and signal simulations to account for alignment and resolution effects that are not fully modeled in the simulations; this procedure will be presented in Section 6. Section 7 describes how we perform a background estimate based on control samples in data, using the lepton |d 0 | to separate signal-like events from background-like ones. Because the lepton |d 0 | distributions for the signal processes are modeled with simulation, we validate the modeling of the displaced tracking efficiency in data using displaced muons from cosmic ray events. Section 8 describes how, from these studies, we derive systematic uncertainties that are applied to the signal.

Event reconstruction and selection
Events were recorded with different triggers in each of the three channels. Because standard CMS electron triggers are not designed to recognize displaced tracks, we use photon triggers instead to ensure efficiency for finding displaced electrons. In fact, photon HLT paths efficiently select electrons as well, as these triggers do not veto electrons or charged particle tracks. In the eµ channel, the trigger requires at least one muon candidate that is not constrained to the primary pp interaction vertex and without any maximum |d 0 | requirement, and at least one photon candidate. In 2016 data, the muon and photon candidates are both required to have p T > 28 GeV if the muon candidate |d 0 | is greater than 0.01 cm, and p T > 38 GeV otherwise. In 2017 and 2018 data, the muon and photon candidates are both required to have p T > 43 GeV; the p T threshold was increased between the 2016 and 2017 data-taking periods to mitigate the effects of the increased pileup. In the ee channel, the events are required to pass at least one of two HLT paths. The first HLT path simply requires at least two photon candidates with p T > 60 (70) GeV in 2016 (2017 and 2018) data. The second HLT path, which is included to partially recover events with lower p T electrons, requires the highest p T photon candidate to have p T > 30 GeV and the second-highest p T photon candidate to have p T > 18 (22) GeV in 2016 (2017 and 2018) data. The photon candidates passing this second trigger must satisfy requirements based on calorimeter cluster shape, isolation, and the ratio of the hadronic to electromagnetic energy, and the diphoton invariant mass must be >90 GeV. In the µµ channel, the trigger requires at least two muon candidates without any primary vertex constraints and without any maximum requirement on the impact parameter. In 2016 data, the muon candidates are required to have p T > 23 GeV if the muon candidate |d 0 | is greater than 0.01 cm, and p T > 33 GeV otherwise. In 2017 and 2018 data, the muon candidates are required to have p T > 43 GeV. For the masses and lifetimes we consider, the efficiency for signal events to pass these triggers is 20-40%, depending on mass, lifetime, analysis channel, and year.
After requiring that the events pass the triggers described above, we preselect wellreconstructed electrons and muons in each channel. Electrons and muons are reconstructed by associating a track reconstructed in the tracking detectors with either a cluster of energy in the ECAL [54,55] or a track in the muon system [56]. The leptons used in this search are reconstructed with the particle-flow (PF) algorithm [57], which aims to reconstruct and identify each individual particle in an event with an optimized combination of information from the various elements of the CMS detector. The primary pp interaction vertex is taken to be the candidate vertex with the largest value of summed physics-object p 2 T . The physics objects used for this determination are the jets, clustered using the jet finding algorithm [58,59] with the tracks assigned to candidate vertices as inputs, and the associated missing transverse momentum, taken as the negative vector sum of the p T of those jets.
In the eµ channel, we preselect events with at least one PF electron and at least one PF muon, while in the same-flavor channels, we preselect events with at least two PF electrons or muons. To retain sensitivity to signals that could produce leptons with the same charge, we make no requirements on the charge of the leptons. In the eµ channel, we require the electrons to have p T > 42 (45) GeV and the muons to have p T > 40 (45) GeV in 2016 (2017 and 2018). In the ee channel, we require the electrons to have p T > 65 (75) GeV, and in the µµ channel, we require the muons to have p T > 35 (45) GeV. These p T requirements ensure that the trigger efficiencies do not depend on p T . In all three channels, we require the electrons and muons to have |η| < 1.5 in order to remove leptons with poorly measured |d 0 |, which would be predominantly at high |η|. Since the signal leptons in the benchmark models are predominantly central, this requirement has a minimal impact on the signal efficiencies. In addition, electrons in the ECAL transition region between the barrel and endcap detectors are rejected because the electron reconstruction in this region is not optimal. This criterion effectively means that electrons are required to have |η| < 1.44. We also reject electrons and muons in certain regions of the η-φ plane, where φ is the azimuthal angle, because two layers of the pixel tracker were not fully functional during certain data-taking periods, resulting in increased |d 0 | mismeasurements. The rejected regions are 1.0 < η < 1.5, 2.7 < φ ≤ π in 2017 and 0.3 < η < 1.2, 0.4 < φ < 0.8 in 2018. This requirement reduces the relative signal efficiency by <1%, and so we apply it for the entire 2017 and 2018 data-taking periods.
Identification requirements, based on energy deposits in the calorimeters and on hit information in the tracker and muon systems, are imposed on the electrons and muons at the "tight" working point [54][55][56]. Included in these identification requirements is the criterion that at least one of the first two functional pixel layers traversed by the electron must register a hit. The identification requirements for muons include that each muon track must have at least one hit in the pixel detector, at least six tracker layer hits, and segments with hits in two or more muon detector stations.
To ensure that electrons and muons are isolated from other particles, we calculate the scalar p T sum of all other particles around the electron (muon) within a cone of ∆R ≡ √ (∆η) 2 + (∆φ) 2 < 0.3 (0.4), correct this sum for contributions from pileup, and define the relative isolation as the ratio of this sum to the electron (muon) p T . For each lepton, the pileup correction term is calculated as the area of the lepton's isolation cone multiplied by the average energy per unit area in η-φ space in the event. By including objects from any vertex in the isolation sum, we allow for the possibility that the displaced lepton is associated with the wrong primary vertex. We require that the relative isolation is <0.10 for muons, <0.0588 for electrons in 2016, and <0.0287 + 0.506 GeV/p T for electrons in 2017 and 2018.
Besides these object-level selections, we also impose several event-level selections. To remove cosmic ray muons in the µµ and eµ channels, we require that there are zero pairs of muons with cos α < −0.99, where α is the three-dimensional angle between the muons. This selection removes back-to-back muons, which is how cosmic ray muons from the Earth's atmosphere appear in the detector. In addition, we require that the relative time between the leading (largest p T ) two muons is not consistent with the timing of cosmic ray muons. To determine the time associated with each muon, we propagate the signal times as measured in the DTs and CSCs to the beam axis assuming the muons are outgoing. We then determine which muon is above the other based on their φ measurements and find ∆t, the time of the lower muon subtracted from the time of the upper muon. Since cosmic ray muons traverse the detector from top to bottom, the lower muon appears later in time than the upper muon, assuming the muons are outgoing from the beam axis, making ∆t negative for cosmic ray muons. On the other hand, muons from the pp collision have similar times as they are both outgoing from the beam axis, and so they have a ∆t distribution centered at 0 ns. We reject events with ∆t < −20 ns if there are at least eight independent time measurements to determine the time of flight of each muon. We also require that the relevant leptons in each channel are separated by ∆R > 0.2. This loose requirement is sufficient to significantly reduce the contribution of the heavy-flavor background from cascade decays of b or c hadrons. To remove pairs of leptons from material interactions, we reject events where the candidate leptons form a good displaced vertex that overlaps with the tracker material, which is measured in Ref. [60]. The vertices are reconstructed with the Kalman vertex fitter [61,62], and a "good" vertex is one with a position uncertainty of less than Table 1: The cumulative efficiencies for t → b signal events to pass the 2018 inclusive SR selection, for several choices of t mass (columns) and cτ 0 (rows). For each entry, the numerator is the weighted number of events passing the SR selection, and the denominator is the total weighted number of generated signal events. The corrections described in Section 6 are applied. 10 µm and a χ 2 /N dof < 20, where N dof is the number of degrees of freedom of the fit.
The events satisfying the preselection criteria for each channel are further categorized using the |d 0 | of the selected leptons. We define a "prompt CR" by requiring the electrons and muons to have |d 0 | < 50 µm. This region is dominated by SM processes with prompt leptons and is used to check that the background simulation accurately reproduces the behavior of the data. The "inclusive SR" is defined by requiring the electrons and muons have 100 µm < |d 0 | < 10 cm. We do not select leptons that are displaced by more than 10 cm to ensure that the leptons originate within the pixel tracker, since the lepton identification criteria require hits in at least one pixel layer.
To remove overlaps among the three channels, events that pass the eµ inclusive SR selection are rejected from the same-flavor channels. Table 1 shows the cumulative signal efficiency for several choices of t mass and cτ 0 . The efficiencies are similar for each year of data taking. The larger signal efficiency in the eµ channel relative to either of the same-flavor channels occurs primarily because two top squarks decaying with equal probability to an electron, a muon, or a tau lepton will produce an eµ final state twice as often as either of the same-flavor final states.

Corrections to the simulation
Several corrections are applied to the MC simulations in order to account for the known differences between simulation and data. The simulation is corrected so that its distribution of pileup interactions matches that of 2016, 2017, and 2018 data. In addition, the trigger efficiency is measured in simulation and in an independent data sample recorded with missing p T triggers, for the three data-taking years and for each analysis channel separately. The ratio of the trigger efficiency in data and simulation is applied as a "scale factor" to each event in the simulated samples. Averaged over the years, the trigger scale factors for the eµ, ee, and µµ channels are 0.94 ± 0.01, 1.00 ± 0.17, and 0.88 ± 0.01, respectively. While these trigger scale factors are derived with samples dominated by prompt events, we also evaluate consistent scale factors when the leptons in these samples are required to be displaced. Scale factors are also applied as functions of lepton p T and η in order to correct the performance of the lepton identification and isolation algorithms [54][55][56].
Corrections are also applied to each lepton in order to make the simulated lepton |d 0 | distributions match those in data, in prompt CRs. These corrections are derived by fitting the simulated background and the data lepton d 0 distributions with Gaussian functions in each channel's prompt CR, and by then comparing the widths of the Gaussian functions. The width of each Gaussian function is largely set by the lepton d 0 resolution, and the discrepancy between the data and background MC simulation lepton d 0 distributions is largely due to an overly optimistic alignment in the simulation, which creates an unrealistically precise d 0 resolution. Therefore, we define σ 2 data = σ 2 bkg + σ 2 correction , where σ data is the width of the function fit to data, σ bkg is the width in background simulation, and σ correction is the additional piece that is needed to make the background simulation and data agree in each channel's prompt CR. We find σ data and σ bkg from the fits, and compute σ correction . Because the fit results are similar in each channel, we average the σ correction derived in the ee and eµ channels for electrons, and in the µµ and eµ channels for muons. We then smear the simulated d 0 distribution by picking random values within a Gaussian distribution centered at 0 and with a width of the average σ correction , and applying these random values as corrections to each lepton d 0 . The average σ correction is 14.8 ± 0.4 (9.2 ± 0.4) µm for electrons and 7.6 ± 0.1 (8.1 ± 0.1) µm for muons in 2017 (2018). No |d 0 | correction is found to be needed for the simulation with 2016 data-taking conditions, as the simulation for 2016 data already has refined alignment conditions that match the data. The corrections are applied to both background and signal MC simulation.

Background estimation
Most leptons resulting from SM processes originate in prompt particle decays. However, displaced leptons can arise from several different sources. Displaced leptons resulting from cosmic ray muons, material interactions, and long-lived SM hadrons are largely rejected by the analysis selection criteria. The tight isolation criterion is particularly useful in rejecting the vast majority of the heavy-flavor background. Nevertheless, leptons from mismeasurements of prompt tracks or from decays of tau leptons, which have a mean cτ 0 value of 87 µm, and long-lived hadrons such as B and D mesons, which have mean cτ 0 values of 500 and <100 µm, respectively, may still appear in the SRs.

The ABCD method
To estimate the number of background events in the SRs, we employ an "ABCD method" based on control samples in data that account for all three significant background sources: mismeasurements, tau lepton decays, and heavy-flavor decays. First, we categorize the events that pass the preselection criteria into four regions (A, B, C, and D) based on each lepton |d 0 | measurement, namely, |d a 0 | and |d b 0 |, as shown in Fig. 2. For the eµ channel, |d a 0 | is defined as the leading electron |d 0 | and |d b 0 | is defined as the leading muon |d 0 |. For the ee (µµ) channel, |d a 0 | is defined as the leading electron (muon) |d 0 | and |d b 0 | is defined as the subleading electron (muon) |d 0 |. We use the number of background events in regions A, B, and C to estimate the expected background in region D, which is the SR. We start from the assumption that N B /N A = N D /N C so that the number of background events in D is N B N C /N A , where N X is the number of background events in the given region. This assumption is valid if |d a 0 | and |d b 0 | are statistically independent. We identify deviations from this assumption using the closure tests described in Section 7.3 and correct for them using the procedure described in Section 7.4.

Signal regions
We subdivide the inclusive SR to define nonoverlapping SRs in |d 0 |: • SR I: 100 < both |d 0 | < 500 µm • SR II: 100 < |d a 0 | < 500 µm, 500 µm < |d b 0 | < 10 cm • SR III: 500 µm < |d a 0 | < 10 cm, 100 < |d b 0 | < 500 µm • SR IV: 500 µm < both |d 0 | < 10 cm We also split SR I, which has the largest number of expected background events, into two bins. In the eµ and µµ channels, these bins are in the leading muon p T , and in the ee channel, these bins are in the leading electron p T . The bin boundary is at 90 (140) GeV for the 2016 (2017+2018) eµ channel, at 300 (400) GeV for the 2016 (2017+2018) ee channel, and at 100 GeV for all years in the µµ channel. The p T bins are chosen such that the bin with higher p T contains less than one expected background event, which maximizes the sensitivity to short lifetimes.
Dividing the inclusive SR in this way separates the expected contribution of different background sources into individual SRs, and gives loose SRs with some amount of background contamination but high signal efficiency and tight SRs with little background contamination but also smaller signal efficiency. The signal efficiency in each SR depends on the lifetime of the LLP, so dividing the inclusive SR into multiple SRs also helps the search to be sensitive to a wide range of LLP lifetimes. Because they are nonoverlapping, we can use these SRs simultaneously in the signal extraction procedure.
We perform a separate ABCD estimate for each SR. When performing the estimates, we subdivide regions B and C into 100-500 µm and 500 µm-10 cm regions to match the SR definitions, and subdivide region A and the 100-500 µm subregions of B and C in p T to match the binning of SR I.

Closure tests in one-prompt/one-displaced sidebands
The background estimation method starts from the premise that the lepton |d 0 | values are independent. The preselection criteria remove one possible source of dependence between the |d 0 | values by ensuring that leptons that share a common displaced vertex do not contribute meaningfully to the SRs, but dependence between the |d 0 | values may still arise from processes in which the leptons originate from the same type of parent particle. Specifically, we find that Z → ττ events in which each tau lepton decays to an electron or muon lead to a dependence between the lepton |d 0 | values, since each electron or muon is produced in the decay of a longlived tau lepton. In principle, processes that produce pairs of b or c hadrons could introduce dependence between the |d 0 | values through this same mechanism, but the lepton isolation criteria ensure a negligible SR contribution from events in which both leptons are produced in b or c hadron decays. We therefore expect the degree of |d 0 | dependence to increase with the fraction of events from Z → ττ. Studies with simulation show that leptons from tau lepton decays contribute significantly from about 100 to 500 µm and peak around 200 µm, so we expect the |d 0 | dependence to appear in this range and peak accordingly. Thus, the background contribution from leptons from tau lepton decays is confined to SR I. The ability to measure such dependence between the |d 0 | values depends on the quality of the |d 0 | measurement. Because |d 0 | resolution is better for muons than for electrons and is better in the 2017 and 2018 data-taking periods relative to the 2016 data-taking period, we also expect the dependence between the |d 0 | values to be more apparent in 2017 and 2018 and to increase with the number of muons in the final state. We observe this dependence between the |d 0 | values in the closure tests described below and correct for it using the procedure described in Section 7.4.
We perform closure tests in sideband regions that are orthogonal to the SRs, in data and background simulation. These sideband regions have one prompt and one displaced lepton. We first perform these closure tests in the region where the prompt (displaced) lepton has a displacement of 30-100 (100-500) µm. We use the estimated and actual yields in several subregions of each sideband region to estimate the ratio of the actual yield to the estimated yield in SR I. This procedure will be described in more detail in Section 7.4. Table 2 shows the resulting ratios in data, background simulation, and background simulation with the Z → ττ events from the Z+jets samples removed. As expected, the ratios in data are frequently greater than unity, which indicates nonclosure of the ABCD method and positive |d 0 | dependence. The data ratios generally agree with those of the full background simulation, which indicates that the source of nonclosure is modeled in the background simulation. When the Z → ττ events are removed from the simulation, we find that the ratios are consistent with unity. Because the full simulation successfully models the observed nonclosure in data and because removing the Z → ττ events results in closure, we conclude that Z → ττ events are indeed the only meaningful source of |d 0 | dependence. We note that the superior |d 0 | resolution of the upgraded pixel detector frequently results in fewer events in the sideband regions in 2017 and 2018 data than in 2016 data, which can lead to larger uncertainties in the 2017+2018 closure test results shown in Table 2, notably in the µµ channel.
We next perform closure tests in sideband regions where one lepton is prompt (30-100 µm) and the other is even more displaced (500 µm-10 cm). Table 3 shows the ratios of the actual yield to the estimated yield averaged over the two one-prompt/one-displaced sidebands in data, background simulation, and background simulation with the Z → ττ events removed. In each case, the resulting ratios are consistent with unity, and data and background simulation agree. Thus, in contrast to the 100-500 µm region tests, the 500 µm-10 cm region tests show no evidence of |d 0 | dependence. We also note that removing the Z → ττ events from the background simulation does not significantly affect the results, which provides further evidence Table 2: Closure test results in data, background simulation, and background simulation with the Z → ττ events removed in the 100-500 µm region. The extrapolated ratios of the actual yield to the estimated yield (averaged over the two one-prompt/one-displaced sidebands) and their statistical uncertainties are given.
Year and channel Data

Background estimate correction and systematic uncertainty
We now define a procedure to account for the |d 0 | dependence in the ABCD method and assign a systematic uncertainty in the estimate. First, as is done in the closure tests, we divide each one-prompt/one-displaced sideband into two subregions in the displaced lepton |d 0 |: (1) the 100-500 µm region, where we find |d 0 | dependence from tau leptons to be significant, and (2) the 500 µm-10 cm region, where we find |d 0 | dependence from tau leptons to be insignificant. The 100-500 µm (500 µm-10 cm) sideband region is used as a CR for SR I (SRs II-IV). We perform closure tests in data in each sideband subregion and use the ratio of the actual to the estimated number of events as a measure of nonclosure.
From the 500 µm-10 cm region tests, we take the largest deviation of the ratio from unity as a systematic uncertainty in SRs II-IV, and apply no correction. This approach produces systematic uncertainties between about 40 and 200% in the background estimates, whose central values range from 0.003 to 3.6 events.
When the displaced lepton |d 0 | is between 100 and 500 µm, we fit the ratio as a function of the prompt lepton |d 0 | with a straight line, where the slope and y intercept are allowed to vary. We extrapolate the prompt lepton fit to 200 µm (within SR I), which is the value where simulation indicates we should expect the largest contribution from tau lepton decays. The mean lepton |d 0 | of the 100-500 µm bin in the background simulation is also approximately 200 µm. We determine the uncertainty in each extrapolated ratio by simultaneously varying the fit parameters according to their 68% confidence interval, while accounting for correlations between the parameters. We take the average of the two extrapolated ratios (one from each one-prompt/one-displaced sideband) and derive a correction and its systematic uncertainty from this average ratio. If the average is greater than unity, we use the average as a multiplicative correction to the background estimate, and we use the uncertainty in the average as a systematic uncertainty in the background estimate. The uncertainty in the average is obtained through simple uncertainty propagation. In this case, we also vary the 200 µm extrapolation point by ±50 µm, which is the approximate range of the tau lepton contribution as a function of |d 0 |, and apply the variation in the resulting correction as an additional systematic uncertainty in the background estimate. If the average is less than or equal to unity, we set the correction equal to unity and use the uncertainty in the average as a symmetric systematic uncertainty about unity. The size of the correction varies between 1.0 ± 0.6 and 4.2 ± 1.8, depending on channel and year.

Closure tests in SRs
To test the full background estimation procedure, we perform closure tests in background simulation in the four SRs, with all corrections and systematic uncertainties derived from background simulation. In these tests, both leptons are displaced. The results of these closure tests in the SRs are shown in Table 4 with the 2016 and 2017+2018 yields combined in each channel. In this table, the actual and estimated yields are given, as opposed to Tables 2 and 3, which display the ratios. The actual yields are compatible with the estimated yields, which indicates that the correction performs as expected and the systematic uncertainties are sufficient to cover any unforeseen dependency.

Additional studies
In addition, we perform several studies to check for other potential sources of background in the SRs. First, to check the significance of material interactions, we invert the criterion that rejects material interactions. After doing so, we find no events in the SRs in data, across all channels and years. Thus, we conclude that there is no significant background from material interactions after the full selection criteria are applied. Second, to test for the presence of cosmic ray muon events, we invert the cosmic ray muon rejection criteria in the µµ channel and scale the number of events by the efficiency of cosmic ray muon events to survive the cosmic ray muon rejection criteria, which is found to be <0.03% from a dedicated data sample of cosmic ray muon events. With this study, we find a negligible number of cosmic ray muon events in data.
To estimate an upper limit on the amount of heavy-flavor background in the SR and to check that this background is covered within the nominal background prediction, several studies are performed. First, we perform the nominal ABCD method while additionally requiring at least one b-tagged jet, using the medium working point of the combined secondary vertices algorithm (version 2) [63]. In the µµ channel, which has the smallest relative SR contribution from mismeasurements and thus the most sensitivity to heavy-flavor backgrounds, we find that the estimate with at least one b-tagged jet is negligible. In the second study of the heavy-flavor background in the SR, we look at samples in which we invert the isolation criterion for events that pass the µµ preselection, for data and simulated background from muon-enriched QCD multijet events. These samples are dominated by muons from decays of b hadrons, and the QCD multijet simulation describes the data well in the region outside of the Z boson peak in the invariant mass distribution. We perform a naive ABCD estimate with the QCD multijet simulation and find no evidence for |d 0 | dependence, which indicates that the nominal background estimation already accounts for the heavy-flavor background. In the last heavy-flavor background study, we estimate this particular background in the SRs from the ratios of the number of events in each SR to the number of events in the prompt CR, from the QCD multijet simulation in the nonisolated region. We multiply these ratios by a normalization factor obtained from the number of QCD multijet simulated events that pass the nominal µµ preselection. Using this approach, we estimate that the heavy-flavor background is about 2 (20)% of the nominal background estimate in SR I (IV), which is well covered within the nominal prediction uncertainties.
To investigate possible SR contributions from displaced leptonic decays of SM hadrons, we examine 2018 data and QCD multijet simulation in the µµ channel with both the muon isolation and ∆R requirements inverted. When the SR muon displacement criteria are applied, this region is dominated by events with dimuon invariant masses near those of the J/ψ and ψ(2S) mesons. Such events likely result from b hadron decays and are therefore already covered by the studies described in the previous paragraph, but we nevertheless estimate an upper bound on their SR contribution as an independent check. Using data in this inverted-isolation, inverted-∆R region, we find the ratios of the number of events in each SR to the number of events in the prompt CR. We multiply these ratios by a normalization factor found from QCD multijet simulated events and find that the SM hadron contribution is less than 0.2% of the nominal prediction in SR I, which is negligible, and about 17% of the nominal prediction in SR IV, which is well covered by the large systematic uncertainty in this SR.

Systematic uncertainties in the signal efficiency
The systematic uncertainty in the background estimation method is the most significant one in the analysis: varying the nuisance parameter by one standard deviation shifts the best fit signal strength by about 5%. The following paragraphs describe the systematic uncertainties that are applied to the signal efficiency.
The integrated luminosities of the 2016, 2017, and 2018 data-taking periods are individually known with uncertainties in the 1.2-2.5% range [64][65][66], which when combined for the data set used in this analysis results in a total uncertainty of 1.8%, the improvement in precision reflecting the uncorrelated time evolution of some systematic effects.
The simulation of pileup events assumes a total inelastic pp cross section of 69.2 mb, with an associated uncertainty of 5% [67]. The systematic uncertainty arising as a result of the modeling of pileup events is estimated by varying the cross section of the minimum-bias events by 5% when generating the target pileup distributions. The pileup weights are recomputed with these new distributions and applied to the simulated events to obtain the variation in the yields in the inclusive SR. The average uncertainty is <1%. We treat these uncertainties as 100% correlated across the three years of data taking.
The trigger efficiency systematic uncertainty is given by the uncertainty in the measured trigger efficiency scale factors. These uncertainties are about 1% for the eµ and µµ channels and 10-19% for the ee channel. The uncertainty is larger for the ee channel relative to the other two channels because there are fewer events available for the efficiency measurement in this channel. In addition, to cover the systematic variation observed in the muon trigger efficiency in signal simulation over the full |d 0 | range, we assign an additional 20% uncertainty. We treat these uncertainties as 100% correlated across the three data-taking years.
The efficiency to reconstruct displaced, isolated, high-p T muons can be measured using cosmic ray muon events, as they also have these properties. The tracking efficiency of displaced muons is measured using cosmic ray muon events in simulation and data, and this efficiency is also used as a proxy for the tracking efficiency for displaced electrons. We take the difference in the mean efficiency between data and simulation as a systematic uncertainty in the signal yield. This uncertainty is 2-14%, depending on the data-taking year. The 2017 and 2018 systematic uncertainties are treated as fully correlated, while the 2016 uncertainty is treated as uncorrelated with the 2017 and 2018 uncertainties, since the pixel detector was upgraded after the 2016 data taking. The choice of how to correlate the uncertainties does not significantly affect the results.
One selection within the muon identification could have some dependence on |d 0 |, namely, the requirement that the muons have at least one pixel detector hit. We find the efficiency of this criterion in simulated cosmic ray muon events and cosmic ray muon events from data, and we apply the difference in mean efficiency between data and simulation as a systematic uncertainty in the signal yield. The average uncertainty is about 16 (32)% in the eµ (µµ) channel. The 2017 and 2018 systematic uncertainties are treated as fully correlated, while the 2016 uncertainty is treated as uncorrelated with the 2017 and 2018 uncertainties, since the pixel detector was upgraded after the 2016 data taking. Although we apply a similar pixel detector hit requirement on electrons, we do not apply a systematic uncertainty for electrons because it would require a sample of displaced electrons in data, which is difficult to obtain and verify. We checked that adding such a systematic uncertainty would not significantly affect the results.
For the two systematic uncertainties derived with cosmic ray muon events, the largest uncertainty is in 2016, compared with relatively smaller uncertainties in 2017 and 2018. This is in part because the 2016 cosmic data sample is much smaller, meaning that the statistical uncertainties are larger in this year. In addition, the uncertainty is reduced in 2017 and 2018 due to the upgrade of the pixel tracker, which allows for more precise tracking measurements.
To find the systematic uncertainty associated with the corrections to the lepton identification and isolation, we fluctuate the lepton scale factors up and down by their uncertainties and observe the change in the simulated event yields in the inclusive SR. The average uncertainty for electrons is about 3% in the eµ channel and about 7% in the ee channel, while the average uncertainty for muons is <1%. We treat these uncertainties as 100% correlated across the three data-taking years.
To find the systematic uncertainty associated with the corrections to the lepton |d 0 |, we fluctuate the lepton |d 0 | corrections up and down by their uncertainties and observe the change in the simulated event yields in the inclusive SR. We find that this uncertainty is negligible in 2017 and 2018, and there is no |d 0 | correction needed for the 2016 simulation.
The systematic uncertainties in the signal efficiency are summarized in Table 5. 3.1-3.9% -eµ channel, muons 0.04-0.05% 0.06-0.07% 0.05-0.06% -ee channel 2.3-2.5% 6.4-7.9% 6.3-7.7% µµ channel 0.09-0.10% 0.14-0.15% 0.11-0.13% Figure 3 shows the expected number of background events and the observed data, with a representative signal overlaid, in each p T bin and each SR, for each channel. Because the electron |d 0 | values are measured less precisely than those of muons, the background estimates are generally greatest in the ee channel, despite its stricter p T requirements. Furthermore, the 2016 data contributes relatively more background events than the 2017-2018 data because of the improved |d 0 | resolution after the pixel detector upgrade. We also note that the background predictions are lowest in SR IV, especially when there is at least one muon in the final state.

Results
The observed number of events is consistent with the predicted amount of background. All data points are observed to deviate by less than ±2 standard deviations from the expected standard model background, for each analysis channel as well as for the channel combination. In the high-p T SR I bin, which is the most sensitive bin for large top squark masses and small cτ 0 values, particularly cτ 0 1 cm, the eµ channel has the largest signal yield relative to the other two channels. As described above, this is because there are twice as many chances to have one electron and one muon, since the top squarks decay to each lepton flavor with equal probability. In this bin, the ee and µµ channel signal yields are similar. In SR IV, which is the most sensitive bin for large top squark masses and long lifetimes, the µµ channel has the largest signal yield for cτ 0 10 cm, relative to the other two channels. This is because the muon reconstruction and selection efficiency is better than that of electrons, which is particularly true at large |d 0 | values. In this bin, the eµ and µµ channels have similar amounts of background, and the ee channel has the smallest signal yield out of the three channels for cτ 0 10 cm. Therefore, for large top squark masses, the eµ channel is the most sensitive for cτ 0 10 cm and the µµ (ee) channel is the most (least) sensitive for cτ 0 10 cm. For a given cτ 0 and mass, the relative distribution of signal events across SRs is similar for all benchmark signals we consider.
We perform a simultaneous counting experiment in each SR bin for most of the interpretations we consider. However, the ee and µµ channels are fit individually to calculate limits on GMSB models with a e or µ NLSP. We set 95% confidence level (CL) upper limits on the product of the signal production cross section and branching fraction to leptons (σB), using a modified frequentist method [68][69][70][71][72]. This approach uses the profile likelihood ratio determined by pseudo-experiments as the test statistic [71] and the CL s criterion. The expected and observed upper limits are evaluated through the use of pseudodata sets. Potential signal contributions to event counts in the SRs and CRs are taken into consideration, as are correlated statistical uncertainties that arise when CR event counts are used to predict the number of background events in multiple SRs. The systematic uncertainties and their correlations are incorporated in the likelihood as nuisance parameters with log-normal probability density functions. The statistical uncertainties in the signal and background estimates are modeled with gamma functions. By comparing the expected and observed cross section limits to the theoretical cross sections at NLO, mass and cτ 0 exclusion limits are set for each of the models we consider.

Summary
A search has been presented for long-lived particles decaying to displaced leptons in protonproton collisions at √ s = 13 TeV at the LHC. With collision data recorded in 2016, 2017, and 2018, and corresponding to an integrated luminosity of 113-118 fb −1 , no excess above the estimated background has been observed. Exclusion limits have been set at 95% confidence level. Top squarks with masses between 100 and at least 460 GeV have been excluded for 0.01 < cτ 0 < 1000 cm, with a maximum exclusion of 1500 GeV occurring at cτ 0 = 2 cm, where cτ 0 is the proper decay length. These exclusions assume that 100% of the top squarks decay to a lepton and a d or b quark, where the lepton has an equal probability of being an electron, muon, or tau lepton. The following exclusions assume that the superpartners of the left-and righthanded leptons are mass degenerate. Electron superpartners with masses of at least 50 GeV have been excluded for 0.007 < cτ 0 < 70 cm, with a maximum exclusion of 610 GeV occurring at cτ 0 = 0.7 cm. Muon superpartners with masses of at least 50 GeV have been excluded for 0.005 < cτ 0 < 265 cm, with a maximum exclusion of 610 GeV occurring at cτ 0 = 3 cm. Tau lepton superpartners with masses of at least 50 GeV have been excluded for 0.015 < cτ 0 < 20 cm, with a maximum exclusion of 405 GeV occurring at cτ 0 = 2 cm. In the case that electron, muon, and tau lepton superpartners are mass degenerate, lepton superpartners with masses between 50 and at least 270 GeV have been excluded for 0.005 < cτ 0 < 265 cm, with a maximum exclusion of 680 GeV occurring at cτ 0 = 2 cm. For sleptons with cτ 0 < 0.8 cm, these are the most sensitive results published to date. For 0.10 < cτ 0 < 12 cm, branching fractions greater than 0.03% have been excluded for 125 GeV Higgs bosons decaying to two long-lived scalar particles, assuming each has a mass of 30 GeV and decays with equal probability to electrons or muons. For scalar particles with 0.1 < cτ 0 < 1000 cm that decay to any final state, these are the most sensitive results published to date.

Acknowledgments
We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid and other centers for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC, the CMS detector, and the supporting computing infrastructure provided by the follow-