Search for heavy resonances and quantum black holes in e$\mu$, e$\tau$, and $\mu\tau$ final states in proton-proton collisions at $\sqrt{s}$ = 13 TeV

A search is reported for heavy resonances and quantum black holes decaying into e$\mu$, e$\tau$, and $\mu\tau$ final states in proton-proton collision data recorded by the CMS experiment at the CERN LHC during 2016-2018 at $\sqrt{s}$ = 13 TeV, corresponding to an integrated luminosity of 138 fb$^{-1}$. The e$\mu$, e$\tau$, and $\mu\tau$ invariant mass spectra are reconstructed, and no evidence is found for physics beyond the standard model. Upper limits are set at 95% confidence level on the product of the cross section and branching fraction for lepton flavor violating signals. Three benchmark signals are studied: resonant $\tau$ sneutrino production in $R$ parity violating supersymmetric models, heavy Z' gauge bosons with lepton flavor violating decays, and nonresonant quantum black hole production in models with extra spatial dimensions. Resonant $\tau$ sneutrinos are excluded for masses up to 4.2 TeV in the e$\mu$ channel, 3.7 TeV in the e$\tau$ channel, and 3.6 TeV in the $\mu\tau$ channel. A Z' boson with lepton flavor violating couplings is excluded up to a mass of 5.0 TeV in the e$\mu$ channel, up to 4.3 TeV in the e$\tau$ channel, and up to 4.1 TeV in the $\mu\tau$ channel. Quantum black holes in the benchmark model are excluded up to the threshold mass of 5.6 TeV in the e$\mu$ channel, 5.2 TeV in the e$\tau$ channel, and 5.0 TeV in the $\mu\tau$ channel. In addition, model-independent limits are extracted to allow comparisons with other models for the same final states and similar event selection requirements. The results of these searches provide the most stringent limits available from collider experiments for heavy particles that undergo lepton flavor violating decays.


Introduction
Extensions of the standard model (SM) can accommodate heavy particles that undergo lepton flavor violating (LFV) decays, motivating thereby searches for deviations from the SM in the eµ, eτ, and µτ final states. This paper reports a search for such phenomena in these final states in proton-proton (pp) collisions at the CERN LHC, and is designed to be as modelindependent as possible. Additionally, the results are interpreted in terms of the characteristics of the following possible states: a τ sneutrino ( ν τ ), which can be the lightest SUSY particle (LSP) [1,2] in R parity violating (RPV) supersymmetric (SUSY) models [3], a heavy Z gauge boson in LFV models [4], and quantum black holes (QBHs) [5][6][7][8]. This is the first analysis searching for high-mass lepton flavor violating signals using the full Run 2 data set, recorded at √ s = 13 TeV.
In RPV SUSY models, lepton flavor and lepton number are violated at lowest (Born) level in interactions between fermions and their superpartners. For resonant ν τ signals, the trilinear RPV part of the superpotential can be expressed as where i, j, and k are generation indices, L and Q are the SU(2) L doublet superfields of the leptons and quarks, and E and D are the respective SU(2) L singlet superfields of the charged leptons and down-like quarks. For simplicity, we assume that all RPV couplings vanish, except for those that are connected to the production and decay of the ν τ , and we consider a SUSY mass hierarchy with ν τ as the LSP. In this model, the ν τ can be produced resonantly in pp collisions via the λ coupling, and can decay either into dilepton final states via the λ couplings, or into quarks via λ couplings. We consider only the final states with two charged leptons. Also, this analysis considers ν τ that decay promptly and are not long-lived [9], which we define as having a transverse displacement less than 1 mm from the production vertex. As long as the λ coupling is larger than 10 −7 , a ν τ of mass 1 TeV will not be considered to be long-lived [3].
An extension of the SM through the addition of an extra U(1) gauge symmetry provides a massive Z vector boson [4]. In our search, we assume that the Z boson has identical couplings to SM particles as the SM Z boson, but that the Z boson can also decay to the LFV eµ, eτ, and µτ final states, each with an assumed branching fraction of 10%. This value is defined in a benchmark scenario commonly used in these searches [10,11].
Theories that invoke extra spatial dimensions can lower the fundamental Planck scale to the TeV region. Such theories also provide the possibility of producing microscopic QBHs [5,6] at the LHC. In contrast to semiclassical thermal black holes that can decay to high-multiplicity final states, QBHs are nonthermal objects, expected to decay predominantly to pairs of particles. We consider the production of spin-0, colorless, neutral QBHs in a model with LFV [12], in which the cross section for QBH production depends on the threshold mass m th in n additional spatial dimensions. The n = 1 possibility corresponds to the Randall-Sundrum (RS) braneworld model [13], and n > 1 corresponds to the Arkani-Hamed-Dimopoulos-Dvali (ADD) model [14]. While the resonant ν τ and Z signals are expected to generate narrow peaks in the invariant mass spectrum of the lepton pair, the distribution of the QBH signal is characterized by a sharp edge at the threshold of QBH production, followed by a monotonic decrease at larger masses. Feynman diagrams for all these three models are shown in Fig. 1. This paper is structured as follows. The CMS detector is briefly described in Section 2. A description of the collision data and the simulated event samples used in the analysis is given in Section 3. The event reconstruction and selection are described in Section 4 and the estimation of SM backgrounds is discussed in Section 5. Systematic uncertainties are described in Section 6, followed by the results and their statistical interpretation in Section 7. The paper is summarized in Section 8.
Tabulated results are provided in the HEPData record for this analysis [20].

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 (HCAL), each composed of a barrel and two endcap sections. Forward calorimeters extend the pseudorapidity (η) coverage provided by the barrel and endcap detectors. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. Events of interest are selected using a two-tiered trigger system. The first level (L1), 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 [22]. 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 [23].
The single-muon trigger efficiency exceeds 90% over the full η acceptance, and the efficiency to reconstruct and identify muons is greater than 96%. Matching muons to tracks measured in the silicon tracker results in a relative p T resolution, for muons with p T up to 100 GeV, of 1% in the barrel and 3% in the endcaps. The p T resolution in the barrel is better than 7% for muons with p T up to 1 TeV [24].
The single isolated electron trigger efficiency exceeds 80% over the full η acceptance of |η| < 2.5, and the efficiency to reconstruct electrons is greater than 95% for electrons with transverse momentum larger than 20 GeV [25]. The electron momentum is estimated by combining the energy measurement in the ECAL with the momentum measurement in the tracker. The momentum resolution for electrons with p T ≈ 45 GeV from Z → ee decays ranges from 1.7 to 4.5%. It is generally better in the barrel region than in the endcaps, and also depends on the bremsstrahlung energy emitted by the electron as it traverses the material in front of the ECAL [26].
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. [27].

Collision data and simulated events
The data sample used for this analysis was collected during 2016-2018 pp operation at a centerof-mass energy of 13 TeV. After applying data quality requirements, such as requiring minimum fractions of active detector channels, the total integrated luminosity is 138 fb −1 .
Simulated samples of signal and background events are produced with several event generators. The RPV SUSY ν τ , Z , and QBH signal events are generated at leading order (LO) precision, using the CALCHEP 3.6 [28], PYTHIA 8.203 [29], and QBH 3.0 [12] Monte Carlo (MC) generators, respectively. The width of the Z signal is taken as 3% of its mass, similar to that of the Z boson. The RPV ν τ and QBH signals are generated using the CTEQ6L [30] parton distribution functions (PDF), while the Z boson signals are simulated using the NNPDF 3.0 PDF set in 2016 and the 3.1 PDF set [31] in 2017-2018, respectively. The LO RPV SUSY ν τ signal event yield is normalized to a next-to-leading order (NLO) calculation of the production cross section [32]; in this calculation the factorization and renormalization scales are set to the mass of the ν τ . The POWHEG v2.0 event generator [33][34][35][36] is used to simulate both top quark pair (tt) production, the dominant background over most of the dilepton mass region, and single top quark production. Diboson production, a significant background in the high mass region, is simulated at LO using the PYTHIA generator.
The MADGRAPH5 aMC@NLO v2.2.2 generator [37] is used to simulate Drell-Yan+jets production simulated at NLO with the FxFx jet matching and merging [38]. The cross sections used to normalize the contribution of these backgrounds are calculated at next-to-next-to-leading order for WW, ZZ, single top quark, and tt processes, and at NLO precision for WZ and Drell-Yan events. The POWHEG and MADGRAPH5 aMC@NLO generators are interfaced with PYTHIA for parton showering, fragmentation, and decays. The PYTHIA parameters for the underlying event description are set to the CUETP8M1 (CP5) tune [39] in 2016 (2017-2018) simulated samples.
The generated events are processed through a full simulation of the CMS detector, based on GEANT4 [40][41][42]. The simulated events incorporate additional pp interactions within the same or nearby bunch crossings, termed pileup, that are weighted to match the measured distribution of the number of interactions per bunch crossing in data. The simulated event samples are normalized to the integrated luminosity of the data. The products of the total acceptance and efficiency for the three signal models in this analysis are determined through MC simulation. The trigger and object reconstruction efficiencies are corrected to the values measured in data.

Event reconstruction and selection
The global event reconstruction is performed using a particle-flow algorithm [43], which aims to reconstruct and identify each individual particle with an optimized combination of all subdetector information. In this process, the identification of the particle type (photon, electron, muon, and charged or neutral hadron) plays an important role in the determination of the particle's direction and energy.
The primary vertex (PV) is taken to be the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone, as described in Section 9.4.1 of Ref. [44].
To reconstruct an electron candidate, energy deposits in the ECAL are first combined into clusters, assuming that each cluster represents a single particle. The clusters are then combined in a way consistent with bremsstrahlung emission, to produce a single "supercluster", which represents an electron or photon. These superclusters are used to seed tracking algorithms, and if a resulting track is found, it is assigned to the supercluster to form an electron candidate.
To reconstruct a muon candidate, hits are first fitted separately to trajectories in the inner tracker detector and in the outer-muon system. The two trajectories are then combined in a global-muon track hypothesis.
Hadronic τ lepton decays (τ h ) are reconstructed from jets, using the hadrons-plus-strips algorithm [45], which combines 1 or 3 tracks with energy deposits in the calorimeters, to identify the τ lepton decay modes. Neutral pions are reconstructed as strips with dynamic size in η-φ, where φ is the azimuthal angle in radians, primarily from reconstructed photons, but also from reconstructed electrons which originate due to conversion of photons. The strip size varies as a function of the p T of the electron or photon candidate.
To distinguish genuine τ h decays from jets originating from the hadronization of quarks or gluons, and from electrons or muons, the DEEPTAU algorithm is used [46]. Information from all individual reconstructed particles near the τ h axis is combined with properties of the τ h candidate in the event to yield separate discriminants for jet, electron, and muon backgrounds. The probability that a jet is misidentified as τ h by the DEEPTAU algorithm depends on the p T and quark flavor of the jet. In simulated events of W boson production in association with jets, this probability has been estimated to be 0.43% for a τ h identification efficiency of 70%. The probability that an electron (muon) is misidentified as a τ is 2.60 (0.03)% for a τ h identification efficiency of 80 (>99)%.

Selection of events for the eµ final state
For the eµ selection, at least one prompt, isolated electron and at least one prompt, isolated muon are required in the event. This minimal selection facilitates a reinterpretation of the results in terms of models with more complex signal topologies than a single eµ pair.
Events that satisfy single electromagnetic-cluster or single-muon triggers, with respective p T thresholds of 175 and 50 GeV for photons and muons, are selected. In 2018, the p T threshold of the photon trigger was raised to 200 GeV. Electromagnetic energy deposited by an electron in the calorimeter activates the photon trigger, and the photon trigger is therefore as efficient as the corresponding electron trigger, while its weaker isolation requirements yield an event sample that can also be used in sideband analyses to estimate the background to the signal. Moreover, using only the muon trigger yields about 90% signal efficiency in masses above 1 TeV. Using the photon trigger together with the muon trigger ensures that the signal efficiency is very close to 100% in the high mass region.
The electron candidate must pass the dedicated selection criteria developed for high-energy electrons [25], which require the energy deposition in the ECAL to be consistent with that of an electron. The electron candidates are required to have p T > 35 GeV and |η| < 2.5. The energy sum in the HCAL within a cone defined by ∆R = √ (∆η) 2 + (∆φ) 2 = 0.15 centered around the electron candidate must be <5% of its energy after it is corrected for jet activity unrelated to the primary interaction. The electron candidate must have a prompt track in the η-φ plane that has at most one missing hit in the pixel tracker and that is well-matched between the tracker and ECAL. The high-energy electron selection also requires electrons to be isolated. The scalar-p T sum of tracks within a cone of radius ∆R = 0.3 around the candidate's direction, excluding the candidate track, is required to be <5 GeV. Moreover, the p T sum of energy depositions in the calorimeters within the same cone, excluding the ECAL supercluster, is required to be <3% of the p T of the candidate [25]. The efficiency of the electron identification requirements, measured in data, is approximately 90% for the highly energetic electrons that are relevant for this search.
Muon candidates are required to have p T > 53 GeV and |η| < 2.4. The muon candidate must pass the high-p T muon identification criteria, which require that the transverse and longitudinal impact parameters of muon candidates relative to the primary vertex must be <0.2 and <0.5 cm, respectively. The track of the muon candidate must have at least one hit in the pixel detector and hits in at least six silicon-strip layers, and must contain matched segments in at least two muon detector planes. To suppress backgrounds arising from muons within jets, the scalar-p T sum of all other tracks in the tracker within a cone of ∆R = 0.3 around the muon candidate track is required to be <10% of the p T of the muon candidate. The relative uncertainty in the p T of the muon track is required to be <30%. The efficiency of the muon identification and isolation requirements, measured in data, is about 90% for muons with p T > 50 GeV.
To prevent a loss in signal efficiency resulting from the misidentification of the sign of the electron or muon charge at large p T , the electron and muon are not required to have opposite charges. Since highly energetic muons can produce bremsstrahlung in the ECAL along the direction of the inner-muon trajectory, such muons can be misidentified as electrons. An electron candidate is therefore rejected if there is a muon candidate track with p T > 5 GeV within ∆R < 0.1 of the electron candidate track. Only one eµ pair is considered per event. When there is more than one eµ candidate, the pair with the highest invariant mass is selected for the analysis. The statistical interpretation is conducted by comparing the shape of the observed invariant eµ mass distribution with those expected for the signal and background.

Selection of events for the eτ and µτ final states
Events are required to have at least one prompt, isolated light lepton (electron or muon) and at least one prompt, isolated τ h .
Events that satisfy single-electron triggers with thresholds of 27, 35, and 32 GeV in 2016, 2017, and 2018, respectively, are selected for the eτ channel. To recover efficiency in the high mass region, events that satisfy single electromagnetic cluster triggers are also selected for this channel. The electron candidate must pass the same high-energy electron selection used for the eµ channel, but with an increased p T threshold of 50 GeV.
The single-muon triggers used in the eµ channel are also used to collect the data samples in the µτ channel. The muon candidate must pass the same high-p T muon identification criteria used for the eµ channel.
The τ h candidate having the largest transverse momentum must have p T > 50 GeV and |η| < 2.3. To reduce the rate of jets, electrons, and muons misidentified as τ h , the τ h candidate is required to satisfy the tight, loose, and tight working points of the respective DEEPTAU discriminator [46].
The low transverse mass (m T ) region is dominated by misidentified τ h events, which have no genuine neutrinos, so the missing transverse momentum is small. Thus, a requirement of m T > 120 GeV is applied which helps to reject misidentified τ h background, where m T is defined as: where denotes the light lepton, i.e., the electron in the eτ final state or the muon in the µτ final state, p miss T is the missing transverse momentum vector in the event, and ∆φ is the difference in the azimuthal angle between p T and p miss T .
A muon veto is applied in the eτ final state to remove overlap with the µτ and eµ final states. This veto rejects events if they contain any isolated muon with p T > 35 GeV, |η| < 2.4, passing the high-p T muon identification criteria and having tracker-based isolation <0.15. Events with a well-separated electron pair are also rejected. A well-separated electron pair is defined as two electrons with ∆R(e, e) > 0.5, each of which has p T > 10 GeV, |η| < 2.5, and passes a very loose working point (>95% efficiency) of the electron identification criteria.
In the µτ channel, in order to avoid overlap with the eτ and eµ channels, events are vetoed if they contain an electron with p T > 35 GeV that passes the high-energy electron selection criteria. Events with a well-separated muon pair are also rejected. A well-separated muon pair is defined as two muons passing the high-p T identification criteria, with p T > 10 GeV, |η| < 2.4, tracker-based isolation < 0.15, and with ∆R(µ, µ) > 0.2.
If there is more than one eτ or µτ pair in an event, then the pair with highest invariant mass is chosen. The leptons forming the pair are not required to have opposite electric charges. The statistical interpretations are conducted by comparing the shapes of the observed collinear mass distributions with those expected for the signals and backgrounds. For an eτ pair or a µτ pair, we use the collinear mass, which provides an estimate of the mass of the new resonance or QBH based on their observed decay products. It is justified by the observation that, since the mass scale of the signal is orders of magnitude higher than that of the τ lepton, the τ lepton decay products are highly Lorentz boosted in the direction of the τ candidate. The neutrino momentum can be approximated to have the same direction as the other, visible decay products of the τ lepton, so the component of p miss T in the direction of the visible τ lepton decay products is used to estimate the neutrino p T . The variable x vis τ , which is the fraction of energy carried by the visible decay products of the τ, is defined as x vis τ = p vis T /(p vis T + p miss T, coll ), where p miss T, coll is the part of the p miss T that is collinear with the τ h p T . The collinear mass m coll can then be derived from the visible mass m vis of the eτ or µτ system to be m coll = m vis / √ x vis τ , where m vis is the invariant mass of the visible τ decay products.

Background estimation
The SM background in the LFV dilepton search includes several processes that produce a final state with two different-flavor charged leptons. For all channels, the dominant background contributions originate from tt production. Other less significant backgrounds originate from diboson (WW, WZ, and ZZ), Z → + − ( = e, µ, τ), Wγ, and single top quark production processes. All these backgrounds are estimated from MC simulation. For the Z → ττ process, all decay modes of τ lepton have been considered. Multijet and W+jets processes also con-tribute due to the misidentification of jets as leptons. These backgrounds are estimated from data.
For the eµ channel, to determine the contributions from W+jets and multijet processes to the m eµ distribution, a control sample in data is defined using jet-to-electron misidentification rate (F e ). This rate is obtained from data, by selecting electron candidates with relaxed selection requirements, using events collected with a single electromagnetic-cluster trigger. Data sidebands obtained by inverting the selection requirements on either the electron isolation or shower shape variables are used to evaluate the contribution of jets passing the full electron selection in the control sample [47]. The jet-to-electron misidentification rate is then defined as the number of jets passing the full electron selection divided by the number of jet candidates in the sample. The rate is quantified in bins of misidentified electron p T and pseudorapidity. The measured jet-to-electron rate is then used to estimate the W+jets and multijet contributions to eµ using data containing muons that pass the single-muon trigger and the full muon selection, and electron candidates satisfying relaxed selection requirements, but failing the full electron selection. Each event is weighted by the factor F e /(1 − F e ) to determine the overall contribution from the jet backgrounds. A correction is applied to account for the fraction of events in the control sample that have genuine leptons, estimated from simulation. The background originating from jets misidentified as electrons is about 10% of the total background.
Background from jets mimicking muons is estimated in a similar way, where identification criteria are loosened for muons in order to create a control sample enriched with jet candidates that are misidentified as muons. Then the jet-to-muon misidentification rate F µ is defined as the number of jets passing the full muon selection divided by the number of jet candidates in the sample. The rate is quantified in bins of misidentified muon p T and η as well. Events in data that have one well-identified electron and one candidate satisfying loose requirements of muon identification are selected, with the weight factor F µ /(1 − F µ ) applied to estimate the jetto-muon misidentification contribution in the signal region. The background originating from jets misidentified as muons is about 2% of the total background. Events with both electron and muon misidentified from jets are not considered since their contribution is expected to be small.
In the µτ and eτ channels, the most significant background after the tt and WW processes comes from W+jets and multijet processes, where jets are misidentified as τ h candidates. This background is estimated using collision data, in a control sample with the same selection as the signal region, except that the m T requirement is inverted to m T < 120 GeV. From this low-m T control sample, two subsamples are constructed to compute the probability for an accompanying jet to be misidentified as a τ h : • Subsample A requires τ h candidates to fail the tight antijet discriminator working point but pass the loose working point.
• Subsample B requires τ h to pass the tight working point.
A factor F τ is calculated from these samples, defined as the ratio of the number of events in subsamples B and A. This factor is calculated as a function of the p T of the τ h candidate, the ratio of the p T of the τ h candidate to the p T of its parent jet, and the pseudorapidity of the τ h . The number of events expected from misidentified τ h in the signal region is estimated from a control sample fulfilling the same criteria as events in the signal region, except that the τ h candidates pass the looser working point but fail the tight working point of the antijet discriminator. Each event is weighted by the factor F τ to determine its effective contribution. A correction is applied to account for the fraction of events in the control sample that have genuine τ h candidates, estimated from simulation.

Systematic uncertainties
Various effects impact the shape and the normalization of the invariant or collinear mass distributions and lead to systematic uncertainties in the signal and background estimates. The uncertainty in the modeling of the invariant or collinear mass distributions reflects three types of systematic effects.
The first type includes those that affect both the shape and normalization of the mass distributions. For all channels, the dominant uncertainties arise from the leading tt and subleading WW backgrounds. They result in a 30-50% uncertainty in the number of tt and WW events at a dilepton mass scale of 2 TeV. The uncertainty in the WW background is estimated from the envelope of the resummed next-to-next-to-leading logarithmic calculation of the soft-gluon contributions to the cross section at NLO, as presented in Ref. [48], using changes in the renormalization and factorization scales by factors of 2 and 0.5. Similarly, the uncertainty in the tt background is estimated by considering the variations in PDF and factorization scales, as discussed in Ref. [49]. Other contributions include the uncertainty in the muon momentum scale, which is approximately 1-2% for 1 TeV muons and depends on their η and φ [50,51]. The uncertainty in the muon efficiency (0.5% for 1 TeV muons) is considered in the eµ and µτ channels [51], and the uncertainty in the electron efficiency (varies between 1-5%, depending on p T and η of the electrons) is considered in the eµ and eτ channels [25]. In the τ channels, the uncertainties in the τ h identification (5% for 1 TeV τ leptons) and τ h energy scale (1.5-4.0% for p T (τ h ) > 100 GeV, depending on the decay mode) are considered [45]. Uncertainties in the electron p T scale and resolution, the muon p T scale and resolution, and the pileup rate are also taken into account, but they have negligible impact on the total background. The uncertainties in the determination of the trigger efficiencies have a very small impact on the expected event yields. The uncertainty associated with the choice of the PDF in the simulation is evaluated according to the PDF4LHC prescription [52].
The jet energy scale is determined with an uncertainty amounting to a few percent, depending on the jet p T and η, using the p T -balance method. This correction is applied to Z/γ * → ee, Z/γ * → µµ, γ+jets, dijet, and multijet events [53]. The resulting effect on signal and background expectations is evaluated by varying the energies of jets in simulated events within their uncertainties, recalculating all kinematic observables, such as p miss T , and reapplying the event selection criteria. The effects of uncertainties in the energy scale of the unclustered particles and jet energy resolution are evaluated in a similar way. These systematic uncertainties affect the shapes as well as the normalizations of the collinear mass distributions.
Uncertainties of the second type directly influence only the normalization of the mass distribution. A systematic uncertainty of 1.6% in the integrated luminosity [54][55][56] is taken for the backgrounds and signals. Among the uncertainties in the cross sections used for the normalization of various simulated backgrounds, the 5% uncertainty in the tt background is the most important. A systematic uncertainty of 50% is applied to the estimate of the misidentified jet background derived from data in all three channels. For τ h final states, this uncertainty is obtained by using the misidentification probability F τ derived from an independent control sample of Z (→ µµ)+jets events for the estimation of misidentified jet background. It is found, especially at high masses (≥1.5 TeV), that the collinear mass distributions obtained in the signal region using the F τ calculated from the two independent control samples agree within 50% and are consistent within the statistical uncertainties. Therefore, an overall 50% systematic uncertainty is assigned to the estimation of this background.
Uncertainties of the third type are associated with limited sizes of event samples in the MC simulation of background processes [57]. In contrast to other uncertainties, they are not correlated between bins of the invariant mass distribution.
Taking all systematic uncertainties into account, the resulting relative uncertainty in the background increases with mass. This relative increase does not significantly affect the sensitivity at large mass values, where the expected number of events from SM processes becomes negligible. All the relevant uncertainties are also taken into account in the extraction of various signals.
All uncertainties are considered to be correlated across the different data-taking years, with the exception of the τ h object-related uncertainties and the uncertainty in the unclustered energy, which are derived from statistically independent sources. No correlation between the different final states is considered, and the analysis results are presented independently for each of the final states.

Results and their interpretations
The mass distributions in the eµ, eτ, and µτ channels, shown in Fig. 2, do not exhibit significant deviations from the expected SM background. The last bins with data events show minor statistical fluctuations in the eτ and µτ channels, which are consistent with the SM expectations within 2 standard deviations, as shown later in the limit plots. The expected mass distributions of signal and backgrounds are taken from simulation, with the exception of backgrounds with jets misidentified as leptons, which are estimated using collision data, as discussed in previous sections. Upper limits on the products of the production cross section σ and branching fraction B are determined using a Bayesian binned-likelihood method [58,59] with a uniform positive prior probability density for the signal cross section. The nuisance parameters associated with the systematic uncertainties are modeled via log-normal distributions for uncertainties in the normalization. Uncertainties in the shape of the distributions are modeled via "template morphing" techniques [60]. A Markov Chain MC method [61] is used for the integration. All limits presented here are at 95% confidence level (CL).
Model-specific limits for the product of the cross section and the branching fraction were obtained for the RPV SUSY, Z , and QBH signals and are shown in Figs. 3, 4, and 5, respectively. These are full Run 2 results based on 2016-2018 data sets. The observed and expected lower mass limits obtained in all three channels are summarized in Table 1      In the narrow-width approximation, the value of σB scales with the RPV couplings, in all three channels. For example, in the eµ channel, the following approximation holds: Using the narrow-width approximation formula for the RPV signal cross section, the cross section limit is translated into exclusion bounds in the plane of mass and coupling of the parameter space of the RPV SUSY model for fixed values of the λ couplings responsible for the decay of the ν τ . Limit contours in the plane of mass and coupling for several fixed values of the coupling are shown in Fig. 6.
Model-independent cross section limits are also obtained. The model-specific shape analysis assumes a certain signal shape in invariant (eµ channel) or collinear mass (τ channels). However, alternative new physics processes yielding the LFV final states could cause an excess of a different shape. A model-independent cross section limit is determined using a single bin with a lower threshold on invariant (collinear) mass, and no upper threshold. No assumptions on the shape of the signal distribution are made other than that of a flat product of acceptance times efficiency, Aε, as a function of the mass. The excluded cross section model-independent limit is shown in Fig. 7. In order to determine the limit for a specific model from the modelindependent limit described here, the model-dependent part of the efficiency and acceptance needs to be applied. The experimental efficiencies for the signal are already taken into account.
A factor f m that reflects the effect of the threshold mass m min on the signal is determined by counting the events with masses m > m min and dividing the result by the number of MCgenerated events. The reconstruction efficiency is nearly constant over the entire mass range probed here, therefore f m can be evaluated at the generator level. A limit on the product of the cross section and branching fraction (σBAε) excl can be obtained by dividing the excluded cross section of the model-independent limit (σBAε) MI given in Fig. 7 by the calculated fraction f m (m min ): Here, B is the branching fraction of the new particle decaying to the relevant LFV final state. Models with a theoretical cross section (σB) theo larger than (σB) excl can be excluded. The fraction of events f m (m min ) must be determined for the particular model under consideration.

Summary
A search has been conducted for heavy particles that undergo lepton flavor violating decays into eµ, eτ, and µτ final states. The search is based on proton-proton collision data at √ s = 13 TeV recorded during 2016-2018 in the CMS detector at the CERN LHC, corresponding to an integrated luminosity of 138 fb −1 . The data are consistent with expectations from the standard model. Lower limits at 95% confidence level are set on the mass of supersymmetric τ sneutrinos at 4.2 TeV in eµ, 3.7 TeV in eτ, and 3.6 TeV in µτ channels. A Z vector boson with lepton flavor violating couplings is excluded for masses below 5.0, 4.3, and 4.1 TeV in the eµ, eτ, and µτ channels, respectively, assuming a branching fraction of 10%. In the context of the Arkani-Hamed-Dimopoulos-Dvali model with four extra dimensions, values of the threshold mass for quantum black hole production less than 5.6, 5.2, and 5.0 TeV are excluded in the eµ, eτ, and µτ channels, respectively. In addition, model-independent limits are provided allowing the results to be interpreted in other models with the same final states and similar kinematic distributions. Limits in the eτ and µτ final states, as well as model-independent limits, are reported for the first time in the context of a high-mass lepton flavor violation search. These are the first results of a high-mass lepton flavor violation search using the full Run 2 data set, and they are currently the most stringent limits from any collider experiment.

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 following funding agencies: [10] CMS Collaboration, "Search for lepton-flavor violating decays of heavy resonances and quantum black holes to eµ final states in proton-proton collisions at √ s = 13 TeV", JHEP 04 (2018) 073, doi:10.1007/JHEP04(2018)073, arXiv:1802.01122.