Search for long-lived heavy neutral leptons with displaced vertices in proton-proton collisions at $\sqrt{s} =$ 13 TeV

A search for heavy neutral leptons (HNLs), the right-handed Dirac or Majorana neutrinos, is performed in final states with three charged leptons (electrons or muons) using proton-proton collision data collected by the CMS experiment at $\sqrt{s} =$ 13 TeV at the CERN LHC. The data correspond to an integrated luminosity of 138 fb$^{-1}$. The HNLs could be produced through mixing with standard model neutrinos $\nu$. For small values of the HNL mass ($\lt$ 20 GeV) and the square of the HNL-$\nu$ mixing parameter (10$^{-7}$-10$^{-2}$), the decay length of these particles can be large enough so that the secondary vertex of the HNL decay can be resolved with the CMS silicon tracker. The selected final state consists of one lepton emerging from the primary proton-proton collision vertex, and two leptons forming a displaced, secondary vertex. No significant deviations from the standard model expectations are observed, and constraints are obtained on the HNL mass and coupling strength parameters, excluding previously unexplored regions of parameter space in the mass range 1-20 GeV and squared mixing parameter values as low as 10$^{-7}$.


Introduction
The discovery of neutrino oscillations [1][2][3] has provided experimental evidence for nonzero neutrino masses [4]. Constraints from cosmological considerations [5,6] and several direct measurements [7] indicate that the neutrino masses are extraordinarily small compared to those of other standard model (SM) fermions. A possible explanation is the existence of new heavy neutral leptons (HNLs) with right-handed chirality. The existence of these hypothetical particles would give rise to gauge-invariant mass terms for the SM neutrinos by means of a see-saw mechanism [8][9][10][11][12][13][14][15]. The HNLs, referred to as N in diagrams and equations hereafter, are singlets under all SM gauge groups and therefore cannot interact with the SM particles through either the electroweak or the strong interaction. They can, however, be produced through mixing with the SM electron, muon, and tau neutrinos [16][17][18], with the corresponding matrix elements denoted as V Ne , V Nµ , and V Nτ , respectively. Additionally, HNLs may be responsible for other unexplained cosmological phenomena. A stable HNL, for example, would make a viable dark matter candidate [19,20]. Furthermore, if at least three generations of Dirac HNLs or two generations of Majorana HNLs exist, CP violation can occur in the HNL system and contribute to the matter-antimatter asymmetry of the universe [21,22].
Many experiments have searched for HNLs in the mass range from a few keV to about 1 TeV [23][24][25][26][27][28][29][30][31][32]. The present search probes the direct production of HNLs in proton-proton (pp) collisions via the decay of W bosons, where the SM neutrino turns into an HNL, and the HNL subsequently decays into either a W boson and a charged lepton , or a Z boson and a neutrino ν. The resulting electroweak gauge boson can then decay leptonically, leading to a final state with three charged leptons and a neutrino: The Feynman diagrams corresponding to the processes described above are shown in Fig. 1 for the case of a W + decay. The HNL decay width in this final state is generally dominated by the W * -mediated diagrams. If the HNL is of Majorana nature, and (or and ν ) can either have the same chirality ( Fig. 1 left) or opposite chirality ( Fig. 1 right). The former process is lepton number violating (LNV), while the latter is lepton number conserving (LNC). In the case of an HNL decay mediated by a W * boson, an LNV decay ( Fig. 1 upper left) can lead to final states with no opposite-sign same-flavor (OSSF) lepton pairs, such as e + e + µ − and e − e − µ + (referred to as e ± e ± µ ∓ ), or µ + µ + e − and µ − µ − e + (referred to as µ ± µ ± e ∓ ). These final states are characterized by relatively small SM background contributions, thus providing a pristine signal for the HNL search. In contrast, decays mediated by a Z * boson ( Fig. 1 bottom) and LNC decays ( Fig. 1 right) are always accompanied by an OSSF lepton pair, resulting in final states such as e + e − µ + and e − e + µ − (referred to as e ± e ∓ µ ± ), or µ + µ − e + and µ − µ + e − (referred to as µ ± µ ∓ e ± ).
To explain neutrino observations and at the same time provide candidates for cosmological observations, HNL models generally need to include three families of HNLs, each of which in principle couples to leptons in all SM generations [33]. A wide range of HNL masses and mixing parameters can give rise to realistic mass values for the light SM neutrinos and explain the matter-antimatter asymmetry in the universe. For example, Ref. [34] finds that an HNL mass scale of 10 GeV and a range of |V N | 2 between 10 -11 and 10 -5 are consistent with these requirements. A large fraction of this phase space is testable at the CERN LHC. For the present study, "simplified" models are considered where only one HNL family couples to exactly one lepton flavor, i.e., only one of V Ne , V Nµ , or V Nτ is nonzero. As a consequence, and (or ν ) Figure 1: Typical diagrams for the production and decay of an HNL through its mixing with an SM neutrino, leading to final states with three charged leptons and a neutrino. The HNL decay is mediated by either a W * (upper row) or a Z * (lower row) boson. In the diagrams on the left, the HNL is assumed to be a Majorana neutrino, thus and in the W * -mediated diagram can have the same electric charge, and lepton number can be violated. In the diagrams on the right, the HNL decay conserves lepton number and can be either a Majorana or a Dirac particle, and, therefore, and in the W * -mediated diagram have always opposite charge. The present study considers only the case where and (or ν ), and in case of the decay mediated by a Z * boson also , are of the same flavor.
always belong to the same lepton generation. The limitations of interpretations based on such simplified models are discussed in Refs. [35][36][37].
The lifetime τ N of an HNL is strongly dependent on |V N | 2 and on the HNL mass m N , and increases rapidly at small masses and low values of the mixing parameter: τ N ∝ m −5 N V −2 N [33]. As a consequence, the kinematic properties and acceptance of HNLs with masses below about 20 GeV are significantly affected by their long lifetimes, resulting in HNL decay lengths ranging from several centimeters to a few meters, depending on the coupling strength. If the HNL has a long lifetime, its decay products ( ± , ∓ , ν ; or ν , ± , ∓ ) emerge from a secondary vertex (SV) that is spatially displaced and distinguishable from the primary vertex (PV) of the pp interaction. Such scenarios are referred to as "displaced" signatures, in contrast to "prompt" signatures where the SV cannot be distinguished, which is the predominant scenario for short HNL lifetimes. The production rates of HNLs also depend on m N and |V N | 2 [38,39].
Previous searches for HNLs performed at the LHC used pp collision data corresponding to approximately 36 fb −1 [29][30][31]. The ATLAS Collaboration has reported on a search for HNLs in the mass range of 4.5 to 50 GeV using prompt and displaced signatures in events with three charged leptons [31]. The search for displaced HNL decays was performed in channels with three muons, or two muons and one electron, providing sensitivity exclusively to |V Nµ | 2 . The CMS Collaboration has performed HNL searches using prompt-lepton final states with two same-sign leptons and one or two jets [30], in a mass range of 20 GeV to 1.7 TeV, as well as in final states with three charged leptons [29], where the HNL production was probed in a mass range between 1 GeV and 1.2 TeV.
The present study is based on the analysis of data collected in 2016, 2017, and 2018 with the CMS detector. The total integrated luminosity of the analyzed data set is 138 fb −1 . Events are selected with three charged leptons (electrons and muons), where one lepton is required to originate from the primary interaction, and the other two leptons are used to reconstruct a displaced vertex. By optimizing the identification of leptons produced in the decay of longlived HNLs and by using a dedicated reconstruction of the displaced HNL vertex, the probed phase space is extended to lower values of the HNL mass compared to the earlier search for prompt signatures in the three-lepton channel [29]. Events with three same-flavor leptons (eee and µµµ) provide sensitivity to HNL decays mediated through W * and Z * bosons. Events with leptons of different flavors (eeµ and µµe, where the first symbol refers to the flavor of the lepton originating from the PV) provide sensitivity only to W * -mediated HNL decays, but can distinguish between LNC and LNV decays and thus differentiate between Dirac and Majorana signatures. The final results of this analysis will be presented as a function of m N and |V N | 2 , separately for HNL couplings to electrons and muons.
Tabulated results are provided in the HEPData record for this analysis [40].

The CMS detector and event reconstruction
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 measured in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. 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. [41].
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 [42]. The second level, known as the high-level trigger, 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 [43].
Charged-particle tracks are reconstructed over the full η range of the tracker, finding charged particles with transverse momentum (p T ) as low as 0.1 GeV, or produced as far as 60 cm from the beam line [44,45]. The silicon tracker used in 2016 measured charged particles within the range |η| < 2.5. For nonisolated particles of 1 < p T < 10 GeV and |η| < 1.4, the track resolutions were typically 1.5% in p T and 25-90  µm in the transverse (longitudinal) impact parameter [44]. At the start of 2017, a new pixel detector was installed [46]. The upgraded tracker measured particles up to |η| < 3.0 with typical resolutions of 1.5% in p T and 20-75 µm in the transverse impact parameter for nonisolated particles of 1 < p T < 10 GeV [47].
Tracks are used to reconstruct the vertices in each event. The candidate vertex with the largest value of summed physics-object p 2 T is taken to be the primary pp interaction vertex. The physics objects used for this determination are the jets, clustered using the jet finding algorithm [48,49] 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.
The global event reconstruction with the particle-flow (PF) algorithm [50] aims to reconstruct and identify each individual particle in an event, with an optimized combination of all sub-detector information. In this process, the identification of the particle type (photon, electron, muon, charged or neutral hadron) plays an important role in the determination of the particle direction and energy. ECAL energy clusters not linked to the extrapolation of any chargedparticle trajectory to the ECAL are reconstructed as photons. A primary charged-particle track associated with potentially many ECAL energy clusters corresponding to this track extrapolation to the ECAL and to possible bremsstrahlung photons emitted along the way through the tracker material is identified as an electron. Muons are reconstructed from tracks in the central tracker consistent with either a track or several hits in the muon system, and associated with calorimeter deposits compatible with the muon hypothesis. Charged-particle tracks neither identified as electrons, nor as muons are attributed to charged hadrons. Finally, neutral hadrons are reconstructed from HCAL energy clusters not linked to any charged-hadron trajectory, or from a combined ECAL and HCAL energy excess with respect to the expected charged-hadron energy deposit.
The electron momentum is estimated by combining the energy measurement in the ECAL with the momentum measurement in the tracker. Within the geometric acceptance of the tracker of |η| < 2.5, 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 [51].
Muons are measured in the range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers. The single muon trigger efficiency exceeds 90% over the full η range, 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 (3)% in the barrel (endcaps) [52].
The missing transverse momentum vector p miss T is computed as the negative vector sum of the transverse momenta of all the PF candidates in an event, and its magnitude is denoted as p miss T [53].

Event simulation
Events from Monte Carlo (MC) simulations are used to derive predictions for signal and several of the background processes. To reproduce the correct multiplicity of pp interactions per bunch crossing (pileup) observed in data, simulated minimum bias events are superimposed on the hard collisions in simulated signal and background events. Simulated events are reweighted to match the distribution of the number of pileup collisions per event in data using the total inelastic pp cross section from Ref. [54]. All events are processed through a simulation of the CMS detector based on GEANT4 [55].
Signal events are generated using the MADGRAPH5 aMC@NLO v2.4.2 generator [56,57] at leading order (LO) in the strong coupling constant α S . The model used to generate signal events extends the SM particle content by introducing up to three right-handed neutrinos [38,[58][59][60]. The production cross section is scaled to next-to-next-to-LO (NNLO) precision using K-factors derived from the comparison of SM samples of W boson production generated at LO precision, using exactly the same settings as the signal samples, to the cross section calculated at NNLO with the FEWZ 3.1 program [61][62][63][64].
In the simulation of signal events, we consider HNL masses m N in the range 1-20 GeV. The HNLs are assumed to couple exclusively to one of the three SM neutrino families with different mixing probabilities |V N | 2 in the range 10 -8 to 10 -1 , depending on the HNL mass. The values of m N and |V N | 2 not only determine the HNL production cross section, but also its mean lifetime and, consequently, acceptance and reconstruction efficiency. For a fixed value of m N , therefore, a simple cross section rescaling is not sufficient to correctly reproduce the behavior of other HNLs with same mass and different |V N | 2 . To emulate an HNL signal sample with a specific value of |V N | 2 , we thus apply per-event weights to the events from all signal samples with the same mass such that the HNL lifetime distribution (taken before the parton shower and detector simulation) matches the predicted distribution for the chosen |V N | 2 value, and additionally we apply a global weight to match the expected cross section.
The signal MC samples are generated assuming Majorana HNLs with both LNC and LNV decays. Dirac HNL scenarios for any (m N ,|V N | 2 ) values are emulated from the Majorana HNL samples, using the fact that the HNL production cross section is the same but the HNL lifetime is twice as large for a Dirac HNL without the LNV decay channels. To that end, all Majorana HNL events are reweighted to produce the expected Dirac HNL lifetime, and all LNV events are classified as LNC events assuming equal selection efficiencies for LNV and LNC events.
Contributions from background processes are estimated mostly from control samples in data. However, MC event simulation is still used to validate the data-driven methods, as well as to derive predictions for minor background processes contributing to the signal regions. Drell-Yan (DY) dilepton production in association with up to four additional jets is simulated at nextto-LO (NLO) in perturbative quantum chromodynamics (QCD) with MADGRAPH5 aMC@NLO and the FxFx algorithm [65] for matching jets from matrix element calculations and parton showers. Simulation of W boson production in association with jets (W+jets) is carried out at LO precision with the MADGRAPH5 aMC@NLO generator using the MLM matching scheme [66]. Top quark and diboson (WW, WZ, ZZ) processes with additional jets are generated with the POWHEG v2 [67][68][69][70][71][72][73][74] and MADGRAPH5 aMC@NLO generators at NLO precision. The production of three vector bosons (W and Z bosons) is simulated with MAD-GRAPH5 aMC@NLO at NLO precision.
The NNPDF3.0 and 3.1 [75,76] parton distribution functions (PDFs) are used for simulating the hard process. The PYTHIA v8.226 and v8.230 event generator [77] is used to model parton shower effects using the CUETP8M2T4 and CUETP8M1 underlying event tunes [78] for the simulated samples produced for the analysis of the 2016 data set, while the CP5 tune [79] is used in the production of the corresponding samples simulated with the 2017 and 2018 detector conditions.

Selection of leptons and jets
Signal events consistent with the HNL decay topology are expected to contain one lepton originating from the PV (prompt lepton), as well as two (displaced) leptons whose tracks reconstruct an SV. The analysis makes use of single-electron and single-muon triggers to select events. The single-electron triggers require the presence of at least one isolated electron with p T > 27 (2016) or 32 GeV (2017-2018 data-taking periods), while the single-muon triggers select events with at least one isolated muon with p T > 24 GeV. An additional set of nonisolated triggers at reduced rate with looser p T requirements on leptons is used for the estimation of backgrounds from misidentified leptons.
Signal leptons are required to be isolated from any hadronic activity in the event. An isolation variable I rel is computed as the scalar p T sum of charged hadrons originating from the PV, neutral hadrons, and photons, all within a cone of ∆R = √ ∆η 2 + ∆φ 2 < 0.3 around the lepton candidate direction at the vertex, divided by the transverse momentum p T of the lepton candidate: The term ρA eff is used to mitigate the contribution of the pileup to the isolation variable: the average p T flow density ρ is calculated in each event using a "jet area" method [80], and the effective area A eff is the geometric area of the isolation cone multiplied by an η-dependent correction factor that accounts for the residual dependence of the isolation on pileup.
Electron reconstruction is based on the combination of ECAL information and a Gaussian sum filter tracking algorithm that accounts for possible bremsstrahlung from the electron [51]. Electrons are required to have |η| < 2.5 to be within the geometric acceptance of the CMS tracking system. To select signal electrons and reduce the rate of misidentified electrons, we impose identification criteria based on the following properties of the reconstructed electron: the electromagnetic shower shape, calorimetric energy ratios, track-cluster matching, track quality, and track impact parameters with respect to the PV.
Prompt electrons are identified using a multivariate (MVA) discriminant [51] with a selection efficiency of about 90%, with p T > 30 (32) GeV in the 2016 (2017-2018) data set. Additional selection requirements are imposed on the transverse and longitudinal impact parameters with respect to the PV of the associated track, |d xy | < 0.05 cm and |d z | < 0.1 cm. Prompt electrons are required to be isolated with I rel < 0.1, and must be matched to the corresponding trigger objects. Prompt electrons passing this selection will be referred to as "tight prompt electrons" hereafter.
Displaced electrons are identified using the standard set of sequential requirements used in prompt electron identification [51], but without requiring a veto on photon conversions, as well as not including the requirement on the maximum number of missing inner hits, since it significantly affects the identification efficiency for electrons not emerging from the PV. "Tight displaced electrons" are required to have p T > 7 GeV and to originate from a nonprimary vertex by requiring |d xy | > 0.01 cm (no |d z | requirement is imposed). In addition, these electrons must be isolated (I rel < 0.2) and satisfy the standard prompt electron identification criteria. Samples enriched in misidentified electrons are selected by applying the tight displaced electron identification criteria with a looser requirement on the relative isolation (I rel < 2). These electrons are referred to as "loose displaced electrons".
Muons are reconstructed by combining the detector information of the tracker and of the muon spectrometer [52]. The geometric compatibility between these separate measurements is used in the further selection of muons. Muons are required to have |η| < 2.4 to fall inside the geometric acceptance of the muon detector. All muons considered for the analysis are identified by loose criteria on the matching between inner tracker tracks and muon spectrometer measurements. Additional loose selection criteria are formulated based on the impact parameters of the tracks with respect to the PV.
Prompt muons must satisfy identification criteria based on the track quality and the matching of the inner tracker track with the measurements in the muon detectors, have p T > 25 GeV, and be matched to the corresponding trigger objects. Prompt muons must be isolated (I rel < 0.1) and associated with the PV by satisfying |d xy | < 0.05 cm and |d z | < 0.1 cm. The muon track in the tracking detector is required to use hits from more than 80% of the layers it traverses. Prompt muons passing this selection will be referred to as "tight prompt muons".
Displaced muons must satisfy a set of standard identification criteria, removing some of the requirements that may reduce the efficiency for reconstructed muons emerging from an SV. "Tight displaced muons" are required to have p T > 5 GeV, not be associated with the PV (|d xy | > 0.01 cm and |d z | < 10 cm), and satisfy a looser isolation requirement (I rel < 0.2). In contrast to the prompt muon identification, no selection is applied on the fraction of hits in the tracking detector, nor on the quality of the matching between the tracks reconstructed in the tracking detector and in the muon system. "Loose displaced muons" are selected with the same identification criteria as used to define tight displaced muons, but with a looser requirement on the isolation (I rel < 2).
We evaluate the efficiencies to reconstruct displaced leptons from an HNL decay in simulated signal samples for typical model parameters. The average efficiency to reconstruct a loose displaced electron with a vertex displacement of 10 (25) cm in the transverse plane is found to be 20-40 (15)(16)(17)(18)(19)(20)%, while the average efficiency for those loose displaced electrons to also satisfy the tight identification requirements is 60-80 (50-60)%. The average efficiency to reconstruct a loose displaced muon with a vertex displacement of 10 (50) cm in the transverse plane is 85-90 (40-50)%, while the average efficiency for those loose displaced muons to also satisfy the tight identification requirements is higher than 95 (80)%. While the selection criteria for displaced electrons and muons are slightly different, the large difference in the observed efficiencies is primarily due to the better resolution and higher efficiency in the basic reconstruction of muons compared to electrons [51,52].
Particles reconstructed by the PF algorithm [50] are clustered into jets using the anti-k T algorithm [48] with a value of 0.4 for the distance parameter, as implemented in the FASTJET [49] package. Jets are required to pass several quality criteria, designed to remove those jet candidates that are likely to originate from anomalous energy deposits. To account for the effect of pileup, nonuniformities of the detector response, and residual differences between jets in data and in simulation, jet energy corrections are applied [81], and the p miss T is modified accordingly. Jets selected for analysis are required to have p T > 25 GeV and |η| < 2.4.
Jets originating from the decay of b hadrons are identified using the DEEPCSV b tagging algorithm [82]. Jets are considered b tagged, and referred to as "b jets" hereafter, if they satisfy the loose selection on the DEEPCSV tagger, corresponding to an efficiency of about 84% in selecting jets containing b hadron decays, with misidentification probabilities of about 11 and 41% for light and c jets, respectively.

Event selection and search strategy
Signal events are characterized by the presence of a prompt lepton (referred to as 1 ), two displaced leptons ( 2 and 3 ), and a neutrino. They are selected with the identification criteria discussed in Section 4. In the following, "lepton" will refer to an electron or muon, and 2 ( 3 ) will refer to the lepton in the pair with the higher (lower) p T . Events are split into final states in which 1 and at least one of 2 or 3 are electrons (eeX) or muons (µµX), as expected for LNC HNL decays. Events with eeX final states are sensitive to |V Ne | 2 , while µµX events are sensitive to |V Nµ | 2 .
Given the low HNL masses considered in this analysis (m N < 20 GeV), 1 has the typical p T spectrum expected for W boson decays, with a Jacobian peak around 40 GeV, while 2 and 3 have rather soft p T spectra, an invariant mass m( 2 3 ) smaller than m N , and a small opening angle. In the absence of significant hadronic activity, 1 and the HNL are typically separated by a large azimuthal angle. In Fig. 2, we show properties of the generated leptons in simulated event samples (i.e., for leptons after the simulation of additional radiation but before detector simulation), for a selection of HNL masses and couplings. These features, along with the possible displacement of 2 and 3 , can be used to identify the two leptons coming from the HNL decay.
The requirement |d xy | > 0.01 cm is imposed on 2 and 3 , since they are expected to originate from the SV. Note that given the rapid variation of the expected HNL displacement as a function of m N and |V N | 2 , this initial displacement requirement is kept loose, and does not fully resolve possible ambiguities between prompt and displaced reconstructed leptons.
Among all the leptons satisfying the prompt selection, the one with the highest p T is chosen as 1 . Among all the selected displaced leptons, the two leptons of any flavor with the highest dilepton p T and opposite sign are selected as 2 and 3 (e ± e ∓ , e ± µ ∓ , µ ± µ ∓ ). If there are no opposite-sign displaced lepton pairs, the event is rejected. This selection requirement reduces background processes with misidentified leptons, while retaining almost full efficiency for the signal.
The decay vertex of the HNL is reconstructed from the tracks associated with 2 and 3 . This is done by fitting the two tracks to a common point with a Kalman-filter approach [83]. The fitted SV provides an estimate for the position of the production vertex of 2 and 3 .
After selecting the events with one lepton from the PV and two leptons forming an SV, the kinematic properties of the low-mass HNL events are investigated and used to separate the signal from the SM backgrounds. Typical backgrounds arise from SM processes containing either one or two misidentified leptons, e.g., tt, single top quark, DY, and W+jets production, or containing a real photon that converts into a lepton pair-mostly electrons-in the detector material, such as Wγ and Zγ.
Compared to these background processes, 2 and 3 in signal events are expected to have a small opening angle, given the small mass and relatively large momentum of the HNL. The ∆R separation between 2 and 3 , ∆R( 2 , 3 ), discriminates the signal from various background processes, especially background events with one misidentified lepton. In order to retain high signal efficiency, a relatively loose requirement of ∆R( 2 , 3 ) < 1 is applied.
In the absence of energetic jets in signal events, 1 is expected to recoil at a large angle in the transverse plane from the HNL, and thus from 2 and 3 . An additional selection criterion is applied to the azimuthal separations between 1 and each of the displaced leptons 2/3 to be The invariant mass of the three charged leptons, m( 1 2 3 ), is limited by the mass of the onshell W boson. Given the relatively low momentum carried away by the neutrino, m( 1 2 3 ) peaks just below the W boson mass, with a steep falloff above 80 GeV and a larger tail at lower masses. We select events with m( 1 2 3 ) values between 50 and 80 GeV. This requirement proves to be particularly efficient against the Zγ background, where the photon radiated by one of the leptons from the Z boson decay undergoes an asymmetric conversion: one of the leptons receives most of the photon momentum, while the other is too soft to be detected or identified. Events with a b jet with p T > 25 GeV are rejected, in order to substantially reduce the background from processes with top quarks. The transverse momentum of the 2 3 dilepton system, p T ( 2 3 ), is required to be larger than 15 GeV, since the low-p T region is dominated by background processes. The decay products of the HNL are emitted at small opening angles, and thus the vector sum of the 2 and 3 momenta should be aligned with the HNL flight direction. We assume that the vector from the PV to the SV corresponds to the HNL flight direction, and calculate the angle between this direction and the dilepton momentum, θ(SV, 2 3 ). Events are required to satisfy cos θ(SV, 2 3 ) > 0.99.
The quality of the reconstructed SV necessarily correlates with the precision on the track parameters of 2 and 3 , as well as the spatial separation between their trajectories at the intersection point. The SV quality is estimated as a probability p SV based on the χ 2 fit from the kinematic vertex fitter. A loose requirement of p SV > 0.001 is applied to remove events that are clearly not compatible with a reconstructed SV. Additionally, we calculate the significance S(∆ 2D ), which is defined as the ratio of the transverse PV-SV distance ∆ 2D and its uncertainty. From simulated event samples, it is found that the distribution of S(∆ 2D ) in signal events is approximately constant for a wide range of mass and coupling values, while background processes typically result in small values of S(∆ 2D ). Thus, a tight requirement of S(∆ 2D ) > 20 is applied to increase the signal-over-background ratio in the signal selection.
In the kinematic region of interest there are several background contributions arising from SM processes that include distinct resonances in the dilepton invariant mass spectrum. Events with an OSSF ( 2 , 3 ) pair are removed when ∆ 2D < 1.5 cm and the invariant mass m( 2 3 ) is consistent with the mass of the ω, φ, J/ψ, or ψ(2S) mesons. The contribution from these processes to the larger displacement region is negligible. In the case where 1 forms an OSSF pair with one of 2 or 3 , we remove events if the invariant mass of that pair is consistent with either the mass of one of the resonances mentioned above, or with one of the higher-mass resonances Υ(1S), Υ(2S), Υ(3S), and Z.
The requirements introduced above represent the baseline selection applied in the analysis and are summarized in Table 1. The selection efficiency for signal events is strongly dependent on the HNL mass and mixing parameters, as well as on the selected final state. The average product of acceptance and efficiency is approximately 0.1% for m N = 2 GeV and |V N | 2 = 8 × 10 -5 , and close to 2% for m N = 12 GeV and |V N | 2 = 10 -6 .
In regions of large displacement, the sensitivity of the analysis greatly benefits from small contributions from the displaced lepton background and from the absence of prompt lepton background contributions. In contrast, signal events corresponding to smaller values of τ N populate regions of smaller displacement, where both the prompt and displaced lepton contributions are still significant. No explicit selection is applied on the SV displacement in order for this search to be sensitive to the full range of τ N values that we consider. The sensitivity of the present search is optimized in bins of the HNL displacement and m( 2 3 ).
3.69 ± 0.08 9.46 ± 0.08 10.02 ± 0.08 10.36 ± 0.08 Z 91.2 ± 10.0 The residual background contributions are dominated by DY and top quark production and processes with photon conversions, and have predominantly small ∆ 2D . In contrast, HNL signal events have typically larger ∆ 2D , depending on m N and |V N | 2 . As a result, ∆ 2D provides the best discrimination between signal and background.
Event categories are defined based on the value of ∆ 2D , with optimized sensitivities for various signal scenarios. Selected events are additionally split into categories corresponding to three different final states per flavor of the prompt lepton. The mixed-flavor final states are split into categories where the prompt and displaced lepton of the same flavor have the same (µ ± µ ± e ∓ and e ± e ± µ ∓ ) or the opposite (µ ± µ ∓ e ± and e ± e ∓ µ ± ) charge, such that the former (latter) categories probe the LNV (LNC) hypothesis. For the same-flavor final states (µµµ and eee, where the two displaced leptons are also always of opposite charge), no distinction between LNC and LNV events is possible, so no additional categorization is performed. The categories are associated with different background rates and process composition, enhancing the sensitivity of the search.
The top quark background processes, where 2 and 3 are typically produced in the semileptonic decay of b hadrons, populate the region m( 2 3 ) < 4 GeV, and are nearly absent at higher dilepton masses. The data are therefore further split into two m( 2 3 ) categories of <4 and >4 GeV. In the high-mass region, only small background contributions are expected originating mainly from DY and W+jets production, and thus fewer ∆ 2D categories are used for m( 2 3 ) > 4 GeV. The event categorization in m( 2 3 ) and ∆ 2D is displayed in Table 2.

Background estimation
The dominant source of background for the final states considered in this search consists of events with misidentified hadrons, muons from pion or kaon decays, and leptons from heavyflavor hadron decays. We refer to such leptons as "background leptons". Events with background leptons originate primarily from top quark, DY, and W+jets production, where the background leptons are typically selected as nonprompt leptons. The contribution of events with displaced background leptons is estimated with a "tight-to-loose" method using control samples in data. In contrast, the contribution of events with prompt background leptons, e.g., from SM events composed uniquely of jets produced through the strong interaction, referred to as QCD multijet events, is found to be negligible, because of the much smaller probability for background leptons to satisfy the tight prompt selection criteria.
Processes that include W or Z bosons can contribute to trilepton final states with at least one displaced electron if a photon radiated from the initial or final state converts into an electronpositron pair in the material of the beam pipe or the detector. This background gives only a minor contribution, and is estimated from control samples in data.
Processes with multiple bosons (W, Z, γ, H) or top quarks can contribute to final states with three or more prompt leptons. The dominant contribution arises from WZ and ZZ diboson production, while processes involving more than two bosons or top quarks have very small production rates and are further suppressed by the lepton and b jet vetoes. The background contribution from all these processes is almost negligible, and is estimated from simulated event samples.
We use data control samples to estimate two classes of background lepton events. Singlebackground (SB) leptons are defined as single reconstructed leptons produced via one of the aforementioned mechanisms. Multiple SB leptons can be found in an event, produced in the decays of different hadrons. In this case, the probabilities of selecting background leptons can be treated as uncorrelated. Double-background (DB) leptons represent pairs of reconstructed leptons produced in the decay chain of the same hadron (e.g., in b → − ν c → − ν + ν s) or a quarkonium state (e.g., in J/ψ → − + ). In such decays, the two leptons are not independent and their selection probabilities are correlated. In background events with background leptons, the main contribution consists of DB leptons, with subdominant contributions arising from SB lepton events.
In the "tight-to-loose" method, the probability f (or "tight-to-loose ratio") for displaced leptons passing a "loose" isolation criterion of I rel < 2 to also pass the "tight" criterion of I rel < 0.2 (which is part of the tight lepton selection described in Section 4) is measured in a "measurement region" selected with requirements orthogonal to those of the baseline event selection. A weight based on this probability is then applied to the number of events found in an "application region", to estimate the background contribution in the signal region. Events in the application region are selected with requirements identical to the baseline selection except for having at least one "loose-not-tight" displaced lepton (i.e., with 0.2 < I rel < 2). This method is applied separately for the estimation of SB and DB lepton events, with separate measurement and application regions in the two cases.
Background leptons are typically found in the proximity of a jet. Each loose displaced lepton is therefore associated with the jet of p T > 10 GeV closest in ∆R, if such a jet exists within ∆R( , jet) < 0.4. If two loose displaced leptons are either associated with the same jet or have ∆R( 2 , 3 ) < 0.45, they are considered DB leptons and treated as a single dilepton system. In all other cases, i.e., if they have ∆R( 2 , 3 ) > 0.45, and are either associated with different jets or at least one is not associated with a jet, the loose displaced leptons are considered SB leptons and treated independently. To reduce the dependency of f on the flavor and momentum of the initial parton from which the background leptons originate, the parton momentum is reconstructed from the background leptons and the associated jets, and f is measured as a function of the reconstructed parton momentum.
For DB leptons, one common parton momentum is reconstructed as the sum of the two displaced lepton momenta and the recalibrated momenta of associated jets. The jet momenta are recalibrated by subtracting the associated lepton momenta from the raw jet momentum, and applying the appropriate jet energy corrections. The tight-to-loose ratio for DB leptons, f DB , is parametrized by the p T of the reconstructed parton. The measurement region for f DB is defined by the baseline event selection summarized in Table 1 except that we also require the presence of at least one b jet and that the two loose displaced leptons pass the DB lepton selection. The b jet requirement enriches the selection in top quark events, which are the main source of DB leptons. We measure the f DB values to be in the range 4-60%, depending on the parton p T . The application region for the DB lepton estimation is also defined by the baseline event selection except for the requirements that at least one displaced lepton is loose-not-tight and that the two displaced leptons pass the DB lepton selection.
For SB leptons, the transverse momentum of the initial parton, p parton T , is reconstructed from the lepton p T and isolation as: where 0.2 is the maximum I rel value allowed for tight displaced leptons. This definition accounts for the hadronic activity in the proximity of the reconstructed lepton, and does not depend on the presence of a jet close to the lepton. The tight-to-loose ratio for SB leptons, f SB , is parametrized as a function of p parton T and of the lepton η. To construct the f SB measurement region we use prescaled single-lepton triggers to select an event sample enriched in multijet events. Events in the measurement region are required to contain exactly one loose displaced lepton, and at least one jet separated from the lepton by ∆R > 1. To reduce the contribution from events with prompt leptons, mainly arising from W+jets and DY production, we require p miss T > 20 GeV and m T = √ 2p miss T p T ( )(1 − cos ∆φ) < 20 GeV, where ∆φ is the azimuthal angle difference between the lepton and p miss T . The residual contribution from prompt lepton events is subtracted using simulated event samples for W+jets, DY, and tt production, with yields normalized to the observed yields in an event selection defined by p miss T > 20 GeV and 70 < m T < 120 GeV. We obtain f SB values of 10-60 (0.1-50)% for electrons (muons), depending on the parton p T and |η|. The application region for the SB lepton estimation is defined by the baseline event selection except for the requirements that at least one displaced lepton is loose-not-tight and that the two displaced leptons do not pass the DB lepton selection.
The accuracy of the background estimates obtained with the SB and DB lepton methods is derived from closure tests performed in data. A background validation region is formed by selecting events that satisfy the baseline event selection summarized in Table 1, except for either fulfilling m( 1 2 3 ) / ∈ [50,80] GeV or having at least one b jet. The contribution of background processes in the background validation region is estimated with the tight-to-loose method, i.e., by applying the measured values of f to events that pass the selection criteria of the background validation region except for having at least one loose-not-tight displaced lepton, and compared to the observed event yields. Using the kinematic regions defined in Table 2, the comparison is shown in Fig. 3. We find excellent agreement within the statistical and systematic uncertainties discussed in Section 8.
Signal contamination in the measurement or application regions used for the SB and DB lepton background contributions would result in an enlarged background prediction and could thus hide the signal contribution in the signal region. We have evaluated this for several HNL signal scenarios that are not excluded by previous searches and find that the signal contribution in the measurement or application regions, when propagated to the signal region in the same way as   and µµX (lower) final states. The hashed band indicates the total systematic and statistical uncertainty in the background prediction. The lower panels indicate the ratio between data and prediction, and missing points indicate that the ratio lies outside the axis range. The uncertainty band assigned to the background prediction includes statistical and systematic contributions. Small contributions from background processes that are estimated from simulation are collectively referred to as "Other".
the background prediction, is negligible compared to the direct HNL contribution in the signal region.

Validation of identification techniques
The contributions from the main background processes in this search are estimated directly from data. The efficiencies for reconstruction, identification, and isolation of displaced leptons in signal events are also validated in dedicated control regions in data.
The performance of the displaced electron reconstruction, identification, and isolation has been studied in data and MC simulation using a Zγ conversion-enhanced selection. The conversion product is taken as a proxy for the displaced lepton. The asymmetric photon conversions include Z → − + γ → − + e ± (e ∓ ) events, where (e ∓ ) represents a low-p T electron that may fail to be reconstructed in the detector. The other three leptons have an invariant mass close to the Z boson mass. Events are selected with two oppositely charged prompt tight electrons with p T > 35 and 7 GeV, respectively, or two oppositely charged prompt tight muons with p T > 28 and 5 GeV, respectively. Additionally, one or two loose electrons are required. The invariant mass of the tight OSSF lepton pair together with one of the additional electrons must be within 10 GeV of the Z boson mass.
While conversion events with four leptons allow for a reconstruction of the SV, and thus directly give a ∆ 2D value, the size of the selected event sample is very small and does not allow for a precise determination of correction factors. The yield of conversion events with three leptons, in contrast, is much larger, but the SV cannot be reconstructed and no ∆ 2D value can be calculated. As a proxy for the actual displacement of the converted photon with respect to the PV, we define the "displacement variable": where r e is the radius of curvature of the path of the displaced electron in the magnetic field, and d xy is the transverse impact parameter of the associated track.
In Fig. 4 (upper left), the distribution of d is shown for conversion events with three leptons. The measured data-to-simulation ratios in bins of d are used as scale factors per electron to apply to the HNL signal events. These scale factors deviate from unity by up to 10% per electron, depending on the displacement value. The four-lepton invariant mass and the ∆ 2D distribution for conversion events with four leptons are shown in Fig. 4 as well, where the scale factors derived from the events with three leptons have been applied to the simulated events. The agreement within uncertainties between the observed and predicted yields validates the correction procedure.
For displaced muons, the efficiencies in the reconstruction, identification, and isolation are studied separately. The general muon reconstruction efficiency is known to be very high and to agree between data and simulation [52], and thus only the displaced-track reconstruction efficiency is evaluated below. Displaced-muon identification efficiencies without the isolation requirement applied are validated using B ± → J/ψK ± → µ + µ − K ± events, where the two muons from the J/ψ meson decay are studied with a "tag-and-probe" method [52]. The measured efficiencies are found to be very close to 100% and the agreement between data and simulation is within the statistical uncertainties. Displaced-muon isolation efficiencies are validated using Z → µ + µ − events, with a tag-and-probe method as well, and no disagreement between data and simulation is observed.  In order to validate the efficiency of the displaced-track reconstruction for muons in simulation, we study K 0 S mesons decaying into two charged pions. To allow for the comparison of the efficiency in data and simulated events, this validation is performed in events with Z → µ + µ − decays. Events are selected with exactly two muons with p T > 30 (25) GeV and with the invariant mass of the muon pairs within 10 GeV of the Z boson mass. Events containing at least one b jet are vetoed to further suppress the top quark background. In these events, K 0 S mesons are reconstructed using pairs of tracks. The number of K 0 S mesons in data and simulation are extracted through a template fit to the invariant mass distribution of the π ± π ∓ pairs, with the combinatorial background subtracted. The results of this validation are presented in bins of radial displacement in Fig. 5.
The event yield in simulation is normalized to the corresponding yield in data at small radial distances (0-0.5 cm), representing an overall normalization correction for possible MC mismodeling of the K 0 S production in Z boson events. The correction scale factors for displaced-track  The invariant mass distribution of the K 0 S candidates reconstructed using the π ± π ∓ tracks for various displacement regions. The fitted K 0 S candidate mass in data in each region is also shown. The lowermost plot shows the K 0 S candidate yield after subtracting the background in data and in simulation, as well as their ratio, as a function of radial distance of the K 0 S vertex to the PV. The simulated yields are scaled to match the data in the lowest displacement bin. reconstruction efficiencies are measured in bins of the radial distance from the PV, as well as in bins of the p T of the reconstructed K 0 S for each data-taking period. In Fig. 5, the comparison of data to simulation for the full data set and all p T ranges is shown as a function of radial displacement. In events with either two displaced muons or one displaced muon, these results are used to extract correction factors for MC signal event acceptances in bins of displacement, as well as the corresponding systematic uncertainties for displaced-track reconstruction. The perevent correction factors deviate from unity by less than 10% for a wide range of displacements and p T . Larger deviations, due to a known inefficiency in the tracker readout electronics [84], are found in the 2016 data set. These only occur for displacements above 4 cm, and result in correction factors deviating from unity by up to 30%.

Systematic uncertainties
Several sources of systematic uncertainties affect the expectations for the number of events in the different signal region bins. The main uncertainties are described below.
The MC generation of HNL signal events is done at LO precision, resulting in large theoretical uncertainties in the cross section of up to 15% when accounting for QCD scale variations and the choice of PDF. To reduce the uncertainty, we apply the K-factor evaluated for W → ν production to HNL signal events. The predicted cross section for the signal process is therefore estimated with NNLO precision, and we obtain a total uncertainty of 3.9% for the combined effects of QCD scale variations and PDF choice.
The integrated luminosity is measured with a precision between 1.2 and 2.5% separately for the three years of data taking [85][86][87]. Because of correlations between some systematic sources in the luminosity calibration, the integrated luminosity of the combined data set is measured with a precision of 1.6%. The uncertainty in the integrated luminosity affects the yields of the signal and of the simulation-based background estimations.
The uncertainty in the modeling of pileup interactions is evaluated by varying the total pp inelastic cross section in simulation by ±5% [54]. This results in an uncertainty in the signal yields of 1-3%, depending on the signal region bin.
To estimate the uncertainties associated with the prompt-lepton identification and isolation efficiency scale factors, the total yields in each signal region bin are recomputed with the scale factors varied up and down by the tag-and-probe fit uncertainties, treating all bins as correlated. Uncertainties of 2-4 (1-3)% in the signal yields are found in events with a prompt electron (muon).
Possible discrepancies in the reconstruction, identification, and isolation efficiencies of displaced leptons between data and simulation are studied using different strategies for electrons and muons, and are used to assess an uncertainty in the MC modeling of the signal. In the different signal region bins, the correction factors applied to the simulated events are varied within their uncertainties, and the resulting uncertainties in the final yields for the displacement bin 1.5 < ∆ 2D < 4 cm are found to be about 7, 4, and 5% for events with ( 2 , 3 ) = ee, eµ, and µµ, respectively.
Trigger scale factors are used to correct the simulated event yields. Uncertainties extracted from the tag-and-probe fits are used to assess systematic uncertainties. These are found to be less than 1% for muon triggers and about 1% for electron triggers.
The uncertainties in the jet energy scale, as well as the jet energy resolution in simulation, affect directly the veto on b jets with p T > 25 GeV. Their impact is estimated by scaling independently the energy of jets up and down by one standard deviation. These systematic effects lead to negligible uncertainties in the analysis.
The efficiency of the b jet veto is corrected for the difference between data and simulation. The uncertainty associated with this correction is used to estimate an uncertainty in the final event yields. The resulting uncertainty is found to be 1-5%.
The contributions from background leptons are estimated using control samples in data, as described in Section 6. The general features of the method were validated in simulation as a function of the search variables. Based on the closure test results in the background validation region in data shown in Fig. 3, we apply a 30% systematic uncertainty in all channels, and include an additional uncorrelated 50% systematic uncertainty in the highest displacement and the m( 2 3 ) > 4 GeV signal region bins, to cover the observed deviations. These two uncertainty sources are considered to be correlated over all data-taking periods and lepton channels. Additionally, we include three uncertainties to account for the different background sources and the differences in statistical power in the f DB measurement regions for the different lepton channels. These uncertainties are in the range 10-17, 10-20, and 25-40%, in the µµ, eµ, and ee channels, respectively, and are taken to be uncorrelated among the channels and years.
The statistical uncertainty in the predicted background lepton contributions in the majority of signal region bins is largely driven by the limited number of events in the application regions. We adopt a Gamma distribution, reflecting the Poisson statistics of the observed event yields in an application region bin with a mean corresponding to the background event yields in the signal region bin divided by f [88,89], to model the statistical uncertainty in the predicted background yields. For each bin, an independent nuisance parameter is assigned accounting for the 68% confidence level (CL) interval of the Gamma distribution with the average measured value of f in that bin. If no events are observed in a given bin in the application region, the Gamma distribution is evaluated with the maximum measured value of f .
For the processes estimated from simulation, the available number of events in the simulated samples limits the precision of the modeling. The statistical uncertainty in the event yield in each search bin is taken as a systematic uncertainty in the shape of the distributions used in the analysis.

Results
The two main variables used to categorize the selected events are ∆ 2D and m( 2 3 ). The comparison between the number of observed events in data and the background predictions is shown in Figs. 6 and 7. The expected and observed yields in each of the signal region bins are shown in Fig. 8, and listed in Tables 3 and 4. In total, 128 (436) data events are observed in the eeX (µµX) final states. The significantly larger event yield in the µµX final states is a consequence of the better resolution and higher efficiency in the reconstruction of muons compared to electrons with the CMS experiment, which correspondingly means that also the signal efficiency is higher in the µµX final states.
The number of observed events in data is in good agreement with the SM background expectations within the statistical and systematic uncertainties and no significant excess is found, for all final states and in all signal region bins. The numbers of expected HNL signal events for a selection of representative signal scenarios are listed in Tables 5 and 6. As discussed in Section 1, only scenarios where the HNL mixes with a single neutrino flavor are considered         Table 5: Number of predicted signal events in the eeX final states. The results are presented for several benchmark signal hypotheses for Majorana HNL production: m N = 2 GeV and |V N | 2 = 0.8 × 10 -4 (HNL2), m N = 6 GeV and |V N | 2 = 1.3 × 10 -6 (HNL6), m N = 12 GeV and |V N | 2 = 1.0 × 10 -6 (HNL12). The quoted uncertainties include statistical and systematic components.  Table 6: Number of predicted signal events in the µµX final states. The results are presented for several benchmark signal hypotheses for Majorana HNL production: m N = 2 GeV and |V N | 2 = 0.8 × 10 -4 (HNL2), m N = 6 GeV and |V N | 2 = 1.3 × 10 -6 (HNL6), m N = 12 GeV and |V N | 2 = 1.0 × 10 -6 (HNL12). The quoted uncertainties include statistical and systematic components. in this search, i.e., these are scenarios where only one of |V Ne | 2 and |V Nµ | 2 is nonzero. The eeX channels (eee, e ± e ∓ µ ± , and e ± e ± µ ∓ ) are sensitive to the |V Ne | 2 mixing parameter, while the µµX channels (µµµ, µ ± µ ∓ e ± , and µ ± µ ± e ∓ ) are sensitive to |V Nµ | 2 .
To derive exclusion limits on HNL signal scenarios, the modified frequentist CL s approach [90][91][92][93] is applied. For each signal scenario, a binned likelihood function L(r, θ) is constructed from the product of Poisson probabilities to obtain the observed yields given the HNL signal scenario scaled with a signal strength r and the background estimates, using the bins defined in Table 2. Additional terms are included to account for all sources of systematic uncertainty, and θ denotes the full set of corresponding nuisance parameters. Based on the profile likelihood method, the LHC test statistic q(r) = −2 ln L(r,θ r )/L(r,θ) is used, whereθ r is the maximum likelihood estimator of θ for a given r, andr andθ are the estimators corresponding to the global maximum of L [92]. The CL s value is calculated in the asymptotic approximation [93], and we exclude the signal scenario if the signal strength value of r = 1 is incompatible with the observed data at 95% CL.
We derive exclusion limits on |V Ne | 2 and |V Nµ | 2 as a function of m N , separately for the cases of Majorana and Dirac HNLs, using a grid of points in the (m N ,|V N | 2 ) parameter space. For m N , we consider values between 1 and 20 GeV, with a step size of 0.5 or 1 GeV. For each mass point, per-event reweighting is used to interpolate to the full range of interest in |V N | 2 starting from the available MC samples. The results are shown in Fig. 9 for a Majorana HNL and in Fig. 10   In the Dirac HNL scenario, we exclude squared values of electron (muon) couplings of 1.1 × 10 -4 (2.0 × 10 -5 ) and higher for an HNL mass of 2 GeV, and between 9.1 × 10 -7 and 1.7 × 10 -4 (3.6 × 10 -7 and 2.3 × 10 -4 ) for an HNL mass of 10 GeV. The highest excluded HNL masses are 14 (16.5) GeV, with |V N | 2 values between 1.1 × 10 -6 and 1.1 × 10 -5 (1.1 × 10 -6 and 2.2 × 10 -6 ). Compared to the previous CMS result targeting prompt signatures [29], the limits are improved by more than one order of magnitude in the probed mass range, which is a result of the dedicated reconstruction of the displaced decay vertex in this analysis. Compared to the results of the DELPHI Collaboration [24], the limits are similarly improved in the probed mass range except for the case of the electron coupling and masses below 3.5 GeV, for which the DELPHI exclusion limits extend to slightly smaller coupling values. Compared to the results of the ATLAS Collaboration [31], which only probe the muon coupling in displaced signatures, we exclude both lower couplings and lower masses. In the Majorana HNL scenario, the exclusion limits are slightly less stringent, with the highest excluded HNL masses being 13 (15) GeV and |V N | 2 values range between 1.1 × 10 -6 and 6.5 × 10 -6 (9.5 × 10 -7 and 1.8 × 10 -6 ) for electron (muon) couplings. The largest explicitly simulated coupling strength is |V N | 2 = 0.5, which is available for the HNL masses of 1, 2, and 3 GeV, and corresponds to a mean lifetime cτ N = 1.3, 0.039, and 0.0045 cm, respectively. We find that this coupling strength is excluded for masses of 1-2 (1-3) GeV in case of electron (muon) couplings.

Summary
A search for heavy neutral leptons (HNLs) has been performed in the decays of W bosons produced in proton-proton collisions at √ s = 13 TeV and collected by the CMS experiment at the LHC. The analysis uses a data set corresponding to an integrated luminosity of 138 fb −1 . Events with three charged leptons are selected, and dedicated methods are applied to identify two displaced leptons consistent with the decay of a long-lived HNL in the mass range 1-20 GeV. Novel methods have been developed to estimate relevant background contributions from control samples in data, addressing one of the important challenges in this type of search at the LHC.
No significant deviation from the standard model predictions is observed. Constraints are derived for models with a single HNL generation of Majorana or Dirac nature, coupled exclusively to electrons or muons. Exclusion limits are evaluated at 95% confidence level on the coupling strengths of HNLs to standard model neutrinos as functions of the HNL mass, covering HNL masses from 1 up to 16.5 GeV and squared mixing parameters as low as 3.2 × 10 -7 , depending on the scenario. These results exceed previous experimental constraints in the mass range 3-14 (1-16.5) GeV for HNLs coupling to electrons (muons) and provide the most stringent limits to date.
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: BMBWF and FWF (Austria); FNRS and FWO (Belgium); CNPq, CAPES, [25] CMS Collaboration, "Search for heavy Majorana neutrinos in µ ± µ ± +jets and e ± e ± +jets events in pp collisions at √ s = 7 TeV", Phys. Lett. B 717 (2012) 109, doi:10.1016/j.physletb.2012.09.012, arXiv:1207.6079.