Heavy Neutrinos with Dynamic Jet Vetoes: Multilepton Searches at $\sqrt{s} = 14,~27,$ and $100$ TeV

Heavy neutrinos $(N)$ remain one of most promising explanations for the origin of neutrinos' tiny masses and large mixing angles. In light of broad advances in understanding and modeling of hadron collisions at large momentum transfers, we revisit the long-standard search strategy for heavy $N$ decaying to multiple charged leptons $(\ell)$, $pp \to N\ell X \to 3\ell \nu X$. For electroweak and TeV-scale $N$, we propose a qualitatively new collider analysis premised on a dynamic jet veto and discriminating, on an event-by-event basis, according to the relative amount of hadronic and leptonic activity. We report that the sensitivity to heavy neutrinos at the CERN Large Hadron Collider (LHC) can be improved by roughly an order of magnitude over the collider's lifetime. At $\sqrt{s}=14$ TeV with $\mathcal{L}=3~{\rm ab}^{-1}$, we find active-sterile mixing as small as $\vert V_{\ell N}\vert^2 = 10^{-2} ~(10^{-3})$ can be probed for heavy Dirac neutrinos masses $m_N \lesssim 1200~(300)$ GeV. The improvement holds also for Majorana $N$ and is largely independent of whether charged lepton flavor is conserved or violated. The analysis, built almost exclusively from inclusive, transverse observables, is designed to be robust across increasing collider energies, and hence serves as a basis for searches at future colliders: With $\mathcal{L}=15~{\rm ab}^{-1}$ at $\sqrt{s}=27$ TeV, one can probe mixing below $\vert V_{\ell N}\vert^2 = 10^{-2} ~(10^{-3})$ for $m_N \lesssim 3500~(700)$ GeV. At a hypothetical 100 TeV $pp$ collider with $\mathcal{L}=30~{\rm ab}^{-1}$, one can probe mixing below $10^{-4}$ for $m_N \lesssim 200$ GeV, below $10^{-3}$ for $m_N \lesssim 4$ TeV, and below $10^{-2}$ for $m_N \lesssim 15$ TeV. We anticipate these results can be further improved with detector-specific tuning and application of multi-variant / machines learning techniques.


Introduction
The origin of tiny neutrino masses and clarity on neutrinos' Majorana nature remain some of the most pressing mysteries in particle physics today. To this end, the hypothesized existence of j ≥ 2 right-handed (RH) neutrinos (ν j R ) [1][2][3][4][5][6][7], represents one of the best, though far from only [7][8][9][10][11][12][13][14][15][16][17], solution to these questions. With RH neutrinos, electroweak (EW)scale Dirac masses can be generated spontaneously during EW symmetry breaking (EWSB), so long as there exists Yukawa couplings to the left-handed (LH) leptons and the Standard Model (SM) Higgs field, L Y = y ij ν L iΦ SM ν j R +H.c. Barring a new fundamental symmetry that ensures lepton number conservation (LNC) below the EW scale, RH neutrinos invariably possess RH Majorana masses µ j R , but possibly also additional Dirac masses (m R ) depending on the precise field content, that can suppress neutrino masses in two distinct limits [18][19][20]: For µ j R much larger than the EW scale v EW = √ 2 Φ SM ≈ 246 GeV, the canonical, highscale Type I Seesaw mechanism [1][2][3][4][5][6][7] is triggered, leading to light neutrino mass eigenvalues inversely proportional to µ j R :M Here, the tilde (∼) denotes matrix objects, withM light representing the light neutrino mass matrix andm D =ỹ ν Φ SM the neutrino Dirac mass matrix. For µ i R that are of the EW scale or much smaller, one realizes the low-scale Type I Seesaw mechanisms, among which are the Inverse [21][22][23][24] or Linear [25,26] Seesaws, and finds light neutrino mass eigenstates with masses proportional to small lepton number violating masses or couplings. For example, in the Inverse Seesaw one hasM In all cases, one predicts the existence of heavy neutrino mass eigenstates (N m ) that can couple to SM particles through mass mixing, perhaps even appreciably. While µ i R is the scale at which lepton number violation (LNV) occurs, without additional inputs, their values in relation to the EW scale are a priori arbitrary. (Though in the vanishing µ R , m D limits, accidental global symmetries are recovered.) Hence, N m may be kinematically accessible at any number of laboratory-based experiments, including collider facilities like the CERN Large Hadron Collider (LHC) and its potential successors [27][28][29][30][31][32][33].
As such, the discovery potential at colliders of particles responsible for neutrino masses, even if only partially accessible, is tremendous [34][35][36][37][38]. More specifically, in the absence of new force carriers and additional sources of spontaneous symmetry breaking, hadron colliders searches for heavy neutrinos cover an impressive range of particle masses, from the MeV up to the TeV, and a multitude of active-sterile mixing angles. (With new gauge bosons and scalars, this can be pushed further.) The limitations of this sensitivity, however, are subtle but substantial: (a) In hadron collisions, projections  and searches [73][74][75][76] closely follow the seminal works of Refs. [34,35,46,77]. Of these, nearly all are based on the resonant production of a single heavy neutrino N m and a charged lepton (ℓ) through the charged Figure 1. Born diagrams for leading heavy neutrino production mechanisms in hadron collisions. Diagrams throughout this text are drawn using JaxoDraw [78].
current (CC) Drell-Yan (DY) process, i.e., qq ′ → W ( * ) → N ℓ, as proposed by Ref. [77]. At the Born level, the process is shown diagrammatically in Fig. 2(a). However, the CC DY process is no longer seen as the only viable means for producing heavy neutrinos at colliders.
It is now known that the W γ → N ℓ ± vector boson fusion (VBF) process [44,55,59,62], shown in Fig. 2(c), is the dominant production mechanism for N m with TeV-scale masses at √ s = 14 TeV [59,62] and sizably enhances the inclusive production rate for lighter N m .
At higher collider energies, neutral current (NC) processes [39,[41][42][43][44], including the gluon fusion (GF) channel [41,42,58] in Fig. 2(b), gg → Z * /h * → N ν ℓ , can surpass the CC DY cross section for masses below 1 TeV [58,66]. For much lighter sterile neutrino masses, the importance of NC and non-resonant processes has likewise been stressed elsewhere [37,38,53,67,71]. (b) The search strategies prescribed by Refs. [34,35,46,77] are highly effective in discriminating against leading SM background processes, but only when reconstructed particles are correctly identified. Misidentified and mis-reconstructed objects, often neglected in theoretical studies, have a substantial impact on sensitivity at hadron colliders [47,59,69,[73][74][75][76]79]. A prime example of this is in multilepton searches for leptonic decays of N m . For production through the CC DY process, the process is given by and shown diagrammatically in Fig. 2. Collider searches for Eq. 1.3 typically require three "analysis-quality" leptons and veto events with additional central, energetic charged leptons. After standard selection cuts, the remaining backgrounds consist of [34,74]: EW processes with many charged leptons, e.g., pp → nW → nℓ + X, wherein one or more charged leptons are too soft or too forward to be readily identified; top quark processes wherein one fails to successfully tag bottom jets; light jets misidentified as electrons or hadronically decaying τ leptons (τ h ); and electrons/positrons whose charges are mis-measured. Undoubtedly, reduction of these so-called irreducible and "fake" backgrounds would have substantial impact on discovery potential.
Born diagram for charged current Drell-Yan production of heavy N with subsequent decay to the three charged lepton (trilepton) final state.
(c) Standard multilepton search strategies for heavy neutrinos are premised almost solely on the existence of high-p T charged leptons originating from heavy N m decays. (This is similarly the case searches of other colorless Seesaw particles [34,38,77].) Aside from vetoing central jets that have additionally been b-tagged, essentially no information about jet activity in exploited in searches for high-mass N m . This is despite the CC DY and W γ fusion mechanisms having qualitatively different QCD radiation patterns than their leading backgrounds. (It is softer, more forward than in leading backgrounds [59,61,69,80].) This is also despite the incredible improvement [81][82][83][84][85][86] in efficiently modeling QCD in hadron collisions for both SM and beyond the SM (BSM) processes since the publication of Refs. [34,35,46,77]. Notably, such improvements include better control over perturbative corrections and uncertainties associated with selection cuts based on the presence of energetic hadronic activity, i.e., jet vetoes [87][88][89][90], in both SM measurements [91][92][93][94][95][96][97] and new physics searches [69,[98][99][100]. Moreover, while never intended for such circumstances, with the incredible increase of top quark and multi-jet cross sections at higher √ s, it is doubtful that current search strategies will remain sufficient at future pp colliders [27,[29][30][31]33]. In light of these observations, we have revisited the long-standard, hadron collider search strategy for heavy neutrinos that decay to fully leptonic final states, as given in Eq. 1.3. We report the construction of a qualitatively new collider analysis that is built almost entirely from inclusive, transverse observables whose shapes remain unchanged with varying collider energies. Unlike past searches, we premise the analysis on the near absence of central, energetic hadronic activity in the CC DY and W γ fusion processes on an event-by-event basis. Using state-of-the-art Monte Carlo (MC) tools, which enable fully differential event generation up to next-to-leading (NLO) in QCD with parton shower (PS) matching, this discrimination is implemented by employing a dynamic jet veto, one where the jet veto p T threshold (p Veto T ) for each event is set to the p T of the leading charged lepton in the event [69]. In doing so, we are able to build an analysis that ultimately selects for the relative amounts of global hadronic and leptonic activities. Importantly, we find irreducible and "fake" background suppression is improved.
For simplicity, we have restricted ourselves to a single, heavy Dirac-like / pseudo-Dirac neutrino as one would find in low-scale Type I Seesaws; for completeness we have also consider a single, heavy Majorana neutrino one would find in baseline phenomenological models. Using as benchmark an analysis inspired by the 2018 CMS experiment's trilepton search for heavy neutrinos [74], we report that the sensitivity to active-sterile mixing (|V ℓN | 2 ) at the 14 TeV LHC can be improved by up to an order of magnitude over the collider's lifetime. By design, the analysis exhibits an unusual robustness against increasing collider energies, and therefore offers to serve as a baseline analysis for future collider searches. We anticipate the results we report can be further improved with detector-specific tuning and application of multi-variant / machines learning techniques [101]. To facilitate such investigations, we make publicly available the MC libraries † used in our computations. The remainder of this report is organized the following order: In Sec. 2 we summarize the key components of the benchmark heavy neutrino models considered in this investigation, commenting also on present constraints. In Sec. 3 is a dedication discussion of our MC modeling, benchmark numerical inputs, and public accessibility of our MC libraries and event samples. Inclusive, parton-and particle-level properties of heavy neutrinos are described in Sec. 4. Jet properties and the theoretical impact of dynamic jet vetoes for both signal and background processes are discussed extensively in Sec. 5. We present the quantitative impact of dynamic jet vetoes in searches for heavy neutrinos at the LHC in Sec. 6. In Sec. 7, we summarize and discuss our findings; we conclude in Sec. 8.

Neutrino Mass Models
If sterile neutrinos participate in the origin of neutrino masses, and if their mixing with active neutrinos is not too small, then it may be possible to produce them at laboratorybased experiments, like the LHC. To model such interactions with SM particles, we consider two benchmark scenarios: The first is the Inverse Seesaw (ISS) model [21][22][23], which is a realistic, low-scale Type I neutrino mass model that permits naturally large Yukawa couplings and heavy neutrinos with masses accessible by colliders experiments. The second is the Phenomenological Type I Seesaw, as investigated by Refs. [34,35]. The qualitative distinction between the two is that the former features pseudo-Dirac / Dirac-like heavy mass eigenstates whereas the latter features heavy Majorana mass eigenstates.
Despite well-established decoupling of EW-and TeV-scale Majorana neutrinos from collider experiments in scenarios where the SM is augmented by only singlet fermions [18][19][20], the latter is considered for two reasons: (i) The discovery of LNV mediated by heavy neutrinos would unambiguously imply a particle field content that is richer than that hypothesized [20]. (This is notable as the canonical collider signatures of the Phenomenological Type I Seesaw can be mimicked by non-minimal scenarios [102][103][104][105][106][107][108]). (ii) More practically, in the absence of new gauge bosons, experimental searches are commonly interpreted in the context of the Phenomenological Type I Seesaw [73][74][75][76].
The remainder of this section is outlined in the following manner: In Sec. 2.1, we discuss the ISS model. In Sec. 2.2, we briefly highlight the relevant distinctions between the Phenomenological Type I Seesaw and the ISS. Constraints on heavy sterile neutrinos are summarized in Sec. 2.3.

Inverse Seesaw
The ISS mechanism [21][22][23] is built on the consequences of the SM being extended by many more SM gauge singlet fermions than the two strictly [109] needed to reproduce oscillation data. The result is a relationship between the scale of LNV and light neutrino masses that is qualitatively different from the traditional high-scale Type I Seesaw. Such a situation can be realized naturally in ultraviolet completions of the SM, including scenarios with warped extra dimensions [110][111][112] and extended gauge symmetries [113][114][115][116][117].
For our purposes, we consider a model wherein the SM field content is extended by three pairs of fermionic singlets, one for each generation of active neutrinos, that possess opposite lepton number, ν j R (L = +1) and X k R (L = −1). A benefit of this field content is that it explicitly displays a nearly conserved lepton number symmetry that is necessarily present [20] in low-scale Seesaw models featuring only singlet fermions. The Lagrangian for the ISS is given by where L SM is the usual SM Lagrangian, L Kin. are kinetic operators for the ν j R and X k R fields, and the so-called "neutrino-portal" couplings and mass terms are given by In Eq. 2.2, Φ SM is SM the Higgs field, and Φ = iσ 2 Φ * is its SU(2) L rotation. Bothỹ ν and m R are complex matrices whileμ X is a complex symmetric matrix whose norm is small compared to the EW scale, in agreement with the near conservation of lepton number. In the ISS, the Majorana masses for ν j R , i.e.,μ R , contribute to light neutrino masses only at a subleading level. They likewise have negligible impact on collider phenomenology. Hence, for simplicity, we takeμ R to be identically zero for the remainder of this study.
After EWSB, the 9 × 9 neutrino mass matrix, in the (ν C L , ν R , X R ) basis, is given bỹ withm D =ỹ ν Φ being the neutrino Dirac mass matrix. Since the mass matrixM ISS is complex and symmetric, it can be put in a diagonal form using the Takagi factorization: . . , m ν 3 , m N 1 , . . . , m N 6 ).

(2.4)
Here m ν 1 , . . . , m ν 3 , are three light neutrino mass eigenvalues, m N 1 , . . . , m N 6 are six heavy neutrino mass eigenvalues, and U ν is a 9 × 9 unitary matrix. Due to the near conservation of lepton number symmetry, the diagonal neutrino mass matrix takes a specific form, with three light neutrino mass eigenstates of Majorana character and six heavy Majorana neutrino mass eigenstates that form three pseudo-Dirac pairs.
For a single generation, i.e., one (ν R , X R ) pair, this can be seen explicitly. In the ISS limit, where µ X ≪ m D , m R , the neutrino mass spectrum takes the form , (2.5) Importantly, one sees the presence of a light mass eigenvalue (m ν ) that is directly proportional to the lepton number violating parameter µ X , which is the inverse of the usual high-scale Type I Seesaw mechanism. One sees also that this LNV scale also controls the mass splitting of the heavier pair of mass eigenstates. Thus, in the limit µ X → 0, i.e., in the limit that lepton number is identically conserved, the light neutrino is massless while the two heavy neutrinos become exactly degenerate, forming a Dirac fermion. Thus, in the context of the ISS, the net contribution of heavy neutrinos N m to lepton number violatingprocesses that occur at a momentum transfer scale Q will be suppressed by the smallness of (µ X /Q) 2 , as required to accommodate light neutrino masses. Moreover, this strongly motivates complementary searches strategies at the LHC (and future colliders) for lepton number conserving processes, such as the trilepton process in Fig. 2 that we consider. With three generations, the neutrino mass matrixM ISS can be diagonalised by block [118], allowing one to write the 3 × 3 light neutrino mass matrix as, which can be further expressed in diagonal form using the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) rotation matrix, U PMNS [119,120]. Doing this gives Importantly, we note that while the PMNS matrix, as historically defined [119,120], is strictly unitary, here, the 3 × 3 block measured by neutrino oscillation experiments is not guaranteed to be separately unitary within the full 9 × 9, unitary mixing matrix and does not strictly correspond to the matrix U PMNS defined above. Indeed, neutrino oscillation experiments measure the sub-block V as defined in Eq.2.9. With these relations, reproduction of oscillation data can be ensured by judicious parametrization of high energy inputs. Specifically, we use the µ X -parameterization [56], giving The above is valid only in the regime that individual elements ofm D are much smaller than individual elements ofm R . For (m ij D /m kl R ) 0.1, including higher order terms in the Seesaw expansion parameter (m Dm −1 R ) is necessary, and can be found in Ref. [121]. In realizing the ISS at collider experiments, we observe first that for the momentum transfer scales Q typically probed in collider experiments, one is insensitive to tiny neutrino masses. In other words, to a very good approximation, (m νm /Q) 2 ∼ 0. This means, however, that on collider scales, the effective Lagrangian one wants should consistently and parametrically follow from the (m νm /Q) 2 ∼ (µ kk X /Q) 2 → 0 limit. For the ISS, this means to zeroth order in (µ kk ′ X /Q) 2 for all k, the relevant heavy neutrino degrees of freedom at collider experiments are the three Dirac mass eigenstates that are recovered precisely in the LNC limit, and henceforth denoted by N m ′ for m = 1, . . . , 3. Now, working in a basis where the charged lepton mass matrix is diagonal, the effective 3 + 3 neutrino mixing matrix U ν in Eq. 2.4 can be generically parameterized by [35] Here, U ℓνm is the 3 × 3 matrix describing light neutrinos (ν m ) coupling to SM weak currents (and measured by oscillation experiments), and V ℓN m ′ is the 3 × 3 matrix describing heavy neutrinos (N m ′ ) coupling to SM weak currents. Subsequently, neutrino flavor states ν ℓ are related to the mass eigenstates by the decomposition [35] where the charged current and charged Goldstone couplings to N m ′ are the analogous neutral current and neutral Goldstone couplings are (2.13) and the coupling to the SM Higgs boson is (2.14) Here, g ≈ 0.65 is the SU(2) coupling constant, θ W is the weak mixing angle, and P L/R = (1 ∓ γ 5 )/2 is the usual LH/RH chiral projection operator for Dirac fermions. With Eq. 2.11, one can determine many properties of heavy Dirac neutrinos, including their partial widths for decays to on-shell EW bosons, which we list here for completeness: The total decay width of N m ′ is then given by the sum over all N m ′ → X partial widths: 18) and likewise the N m ′ → X branching rate (BR) is defined as the ratio .
For fixed |V ℓN m ′ |, in the high-energy / infinite m N m ′ limit, the Goldstone Equivalence Theorem becomes manifest, and one obtains the branching rate relationship [35] BR Contrary to common belief, the 2 : 1 : 1 ratio holds even for Dirac neutrinos due to the fact that the N m ′ − W charged current coupling is a factor of √ 2 larger than the neutral current and Higgs couplings at the Lagrangian level, as seen above. This itself is a manifestation of the fact that the W µ boson is a linear combination of the W 1 µ and W 2 µ vector fields.

Phenomenological Type I Seesaw
The formalism of the Phenomenological Type I Seesaw is well-documented [35] and will not be repeated here in the interest of brevity. For three heavy Majorana neutrinos, the relevant interaction Lagrangian almost identically resembles Eq. 2.11. The crucial difference is the Majorana condition, which states that in the mass basis the field N m ′ and its conjugate are related by N c m ′ = N m ′ . This of course means that, if sufficiently heavy, both the N m ′ → W + ℓ − , Zν ℓ , hν ℓ and N m ′ → W − ℓ + , Zν ℓ , hν ℓ decay channels are open. By CPinvariance, the partial widths for N m ′ → W − ℓ + and N m ′ → W + ℓ − are equal, as are the analogous Z and h partial widths. Hence, the total width for a heavy Majorana neutrino and a heavy Dirac neutrino of the same mass are related by (2.21) This preserves the well-known branching rate relationship,

Constraints on Heavy Sterile Neutrinos
Heavy neutrinos can manifest themselves in a variety of terrestrial-and space-base experiments. Hence their existence is constrained from a number of sources. For reviews on sterile neutrino constraints, see Ref. [38,48,[122][123][124]. We now summarize such constraints: • Global Constraints on Non-Unitarity of PMNS Matrix: Deviation from unitarity of the light neutrino mixing matrix U induced by active-sterile mixing can be expressed generically by the Hermitian matrix η ℓℓ ′ and the relationship In a general Seesaw scenario, constraints from a global fit [124] to EW precision data, tests of CKM unitarity, and tests of lepton universality limit η ℓℓ ′ to be: at 95% CL. In a general Seesaw scenario, η ℓℓ ′ is related to mixing matrices V ℓN m ′ by the relationship Under the benchmark assumption that ℓℓ ′ mixing is dominated/saturated by the lightest heavy mass eigenstate N m ′ =4 , as we invoke, the relationship simplifies to 2|η ℓℓ ′ | ≈ |V ℓ4 V * ℓ ′ 4 |.
• Perturbativity Bounds on Total Width: Large Yukawa couplings, in the form of large active-sterile mixing, can lead to too large heavy neutrino total decay widths.
In this analysis, we will constrain the total width of N m ′ to be:  [59,62] and references therein ‡ . † A technical exception is a change in the heavy neutrino MC particle identification (PID) codes. These differ to avoid conflict with standard HEP PID convention established in 2002 for Majorana neutrinos [132].
‡ For both sets of libraries, model variants that include six massive quarks, three massive charged leptons, and/or non-diagonal quark mixing are also available.
Parton-level event generation at LO and NLO in QCD is handled by mg5amc. For decays of heavy neutrinos and SM particles W, Z, t, we invoke the narrow width approximation (NWA) and pass partonic events to MadSpin [136]. Spin correlation in decays of intermediate resonances are preserves through the implementation of the spin-correlated NWA in MadSpin. Parton events are inputed to Pythia v8.230 (PY8) [137] for parton showering and hadronization. Decays of B hadrons and τ leptons to lighter hadrons and leptons are channeled through PY8. Following the CMS analysis of Ref. [138], we use the PY8 tune CUETP8M1 "Monash*" (Tune:pp = 18) [139]. In addition, for realistic shower simulation, QED shower flags and recoil/primordial momentum flags are switched to on. The impact of underlying event is neglected and left to follow-up investigations. Hadronlevel events are given to MadAnalysis5 v1.6.33 [140,141] for particle-level clustering via FastJet v3.2.1 [142]. Following the jet veto study of Ref. [100], we employ the anti-k T algorithm [143] with radius parameter R = 1. During jet clustering, tagging efficiency for b-jets and hadronic decays of τ leptons (τ h ) is taken to be 100% with 0% misidentification rates of light jets. (Realistic (mis)tagging efficiencies and (mis)identification rates are applied at the analysis-level.) Reconstructed / clustered, particle-level events are ultimately written to file in the standard Les Houches Event (LHE) format [144]. Only jets with p T > 1 GeV are recorded. Particle-level events are analyzed using in-house, ROOT-based code.
In several instances, we compute cross section normalizations that are beyond the accuracy of NLO in QCD with PS/LL matching. For such cases, we heavily exploit the factorization properties within the formalism of Soft-Collinear Effective Field Theory (SCET) [145][146][147], particularly when working in momentum space [148]. Approximate next-to-next-tonext-leading logarithmic (N 3 LL) threshold corrections to the the gluon fusion process in Fig. 1(b), which captures the dominant contributions up to next-to-next-to-leading order (N 2 LO) in α s [149], are obtained following Ref. [66,148,150,151]. In addition, we obtained cross sections with a static jet veto up to NLO in QCD with jet veto resummation matching at N 2 LL, using the automated MadGraph5_aMC@NLO+SCET libraries developed by Ref. [85,96].

Heavy Neutrino Signal Event Generation
We now describe our signal process event generation using mg5amc. Here, we also provide instructions for using the HeavyN_Dirac and HeavyN libraries.
In this study we consider heavy Dirac and Majorana neutrinos produced through the DY and W γ fusion mechanisms, as shown in Fig. 1(a) and (c), and decayed to fully leptonic final states. At the Born level, this corresponds to the production and decay chains, For Dirac neutrinos, the syntax of the HeavyN_Dirac model file is designed to closely match the original HeavyN model file for Majorana neutrinos. Hence, much of the syntax presented in Ref. [62] to simulate Eqs. 3.1-3.2 for Majorana neutrinos is unchanged. For the inclusive CC DY process, the appropriate mg5amc syntax is import model SM_HeavyN_Dirac_NLO set gauge Feynman define p = 21 1 2 3 4 -1 -2 -3 -4 define ww = w+ wdefine ell+ = e+ mu+ ta+ define ell-= e-mu-tadefine ell = ell+ elldefine nn = n1 n1g enerate p p > nn ell QCD=0 [QCD] output HeavyNDirac_DYX_Nl_NLO ; launch -f Here, p here denotes an active parton p ∈ {q, q, g}, with q ∈ {u, c, d, s}, and ell is any one of the three charged leptons, ℓ ∈ {e, µ, τ }. Due to the large top quark mass, neither b or t are considered active partons. The [QCD] flag activates the MC@NLO [152] and MadLoop [153,154] formalisms in mg5amc, allowing one to compute processes up to one loop in α s and parton shower (PS) matching to LL accuracy [85]. For the VBF channel, the corresponding commands are generate p a > nn ell p QED=3 QCD=0 [QCD] add process a p > nn ell p QED=3 QCD=0 [QCD] To make explicit, we consider collinear, initial-state photons from both hadronic and partonic sources, as advised by Refs. [59,62,155,156] and as implemented in Refs. [62]. This is handled automatically by our use of the LUXqed formalism [157,158]; see Sec. 3.4. For the production of heavy Majorana neutrinos, we use the original HeavyN libraries. Two notable (though hopefully obvious) syntax modifications to those above are needed: (i) One needs import the HeavyN model file instead of the HeavyN_Dirac libraries. (ii) One must omit the charge conjugate of N since they are the same particle. Following Ref. [ As the signal process proceeds through an s-channel heavy neutrino, which may be off-shell, we use the commands The operator syntax > X > selects only for diagrams/amplitudes with particle(s) X appearing in the s-channel; || acts as the Boolean inclusive OR operator. The decay of the W boson is treated with the NWA. In both model files, the mass eigenstates N m couple to all SM lepton flavors (ℓ) with a mixing strength of |V ℓNm |. Hence, the above production-level commands will generate matrix elements (ignoring charge multiplicity) for all three (N ℓ 1 )-flavor permutations. Likewise, the decay-level syntax will generate matrix elements for all (ℓ 1 − ℓ 2 − ℓ 3 )-flavor permutations. To help ensure consistent event generation across the various flavor permutations considered, including to correctly account for leptonic decays of τ leptons, we choose to control the relative abundance of various charged lepton flavors by tuning the parameters |V ℓ 1 N 1 |, |V ℓ 2 N 1 |, |V ℓ 3 N 1 | at runtime.
For fully inclusive DY signal samples, no generator-level cuts are needed or are applied. For VBF signal samples, t-channel QED poles involving both quarks and charged leptons are present [59]. These divergences are regulated with the loose generator-level cuts p j, gen. T > 20 GeV, p ℓ, gen. T > 10 GeV, |η j, gen. | < 5, |η ℓ, gen. | < 3. (3.4) As we are considering pp collisions up to 100 TeV, a brief comment is due: For the heavy N mass scales under consideration (m N < 20 TeV), the p Gen.
A final technical comment: mg5amc is less-than-optimized for phase space integration of VBF topologies, particularly at NLO in QCD. This is partially due to a virtual pentagon diagram that, while nonzero in its own right, does not contribute due to color conservation and the lack of color flow at the Born level [163]. To aid computation efficiency, the quantity #IRPoleCheckThreshold in the file mg5amc runtime card Cards/FKS_params.dat is set to -1.0d0 from its original value of 1.0d-5 for all VBF signal computations. For the DY signal process at √ s = 100 TeV, a value of 1.0d-3 or 1.0d-2 is used for m N = 150 − 300 GeV to ease numerical stability requirements for so low Bjorken-x.

Standard Model Background Event Generation
We now elaborate on the details of our SM background event simulation.
To study the trilepton process pp → 3ℓX in conjunction with a jet veto, we consider three-and four-charged lepton processes from a variety of SM sources. In order to better ensure a reliable description of the leading jet at all p T scales, all processes are considered at NLO+PS, with one exception. We separate the SM processes into four categories: (i) top quarks, (ii) EW diboson continuum, (iii) resonant EW multi-boson production, and (iv) fake charged leptons. Throughout the text we refer to "continuum" processes as those that possess resonant and non-resonant components. For such samples, full interference between resonant and non-resonant regions of phase space is considered.

Top Quarks
Due to their large rates and inherent mass scales, top quarks produced in association with EW bosons, e.g., ttW , are leading backgrounds to traditional trilepton searches that rely only on vetoing central, b-tagged jets [34,74]. Invariably, the rates of these many-lepton processes are able to compensate for the small (but not entirely rare) likelihood that b-jets fail to be identified while simultaneously one or more charged leptons fail to meet "analysis quality object" criteria. To account also for these potentially soft (low-p T ), forward (high-|η|) charged leptons, we consider the continuum processes As for the signal process above, p here denotes a parton p ∈ {q, q, g}, with q ∈ {u, c, d, s}, and ℓ ∈ {e, µ, τ }. We consider only pure gauge interactions and neglect the contribution of a SM-like Higgs, e.g., pp → tth, h → τ + τ − . Inclusion of such processes would incredibly complicate our matrix element computations at NLO while contribute only marginally to our background modeling. Owing to the fact that the ttℓℓ → 2ℓX and ttℓℓ → 3ℓX signals are efficiently cut out by minimal selection cuts, for MC efficiency, only the fully leptonic ttℓℓ → 4ℓX decay mode is considered. For the first two top quark processes, which we model at NLO+PS, the mg5amc syntax is: For both processes, MadSpin syntax is: decay tt > ww bb, ww > ell vv The single top process bq → tq ′ ℓℓ is a somewhat novel channel as it was only recently observed at Run II of the LHC [164,165]. Modeling this process is nuanced as we are formally working a four-flavor scheme with a massive b quark (in order to preserve modeling of τ lepton and B meson decays). Hence, for tqℓℓ, we take m b = 10 −10 GeV at the matrix element level in order to use the b quark PDF and circumvent internal mg5amc checks for massive, initial-state partons. (These checks are necessary to ensure consistency of the Collinear Factorization Theorem [166][167][168][169][170][171] as implemented in mg5amc.) On the other hand, standard implementation of the aMC@NLO formalism in mg5amc requires that the initialstate parton b be identically massless, not approximately massless. Subsequently, this is the only event simulation that we model at LO+PS (all other are simulated at NLO+PS). Subsequently, the Born-level tqℓℓ process is modeled with the following: define qq = u c d s u~c~d~sd efine bb = b bd efine tt = t tg enerate qq bb > ell ell qq tt / h QCD=0 QED=4 add process bb qq > ell ell qq tt / h QCD=0 QED=4 The leptonic decay of t is handled in the same manner as the pp → ttℓX processes.
With massive top quarks, the pp → ttℓν matrix element possesses no IR poles at the Born level. Hence, no regulating cuts are needed nor are any applied during event generation. The pp → ttℓℓ and bq → tq ′ ℓℓ matrix elements, on the other hand, do possess IR poles as m ℓℓ → 0. Hence, the following generator-level cuts are applied: m(ℓ 1 , gen., ℓ 2 , gen.) > 10 GeV and |η ℓ, gen. | < 3. (3.8) The invariant mass boundary is set to 10 GeV in order to enrich the contribution of softleptons but high enough to remove contribution of low-mass DY resonances.

Electroweak (Diboson) Continuum
Electroweak production of three-and four-charged lepton events in hadron collisions, pp → ℓℓℓν and pp → ℓℓℓℓ, and zero-resonant (V * V ′ * /V * γ * ) diboson production, are the largest, pure-EW background to trilepton searches for heavy neutrinos with m N at or above the EW scale. In the presence of a jet veto, the 3ℓν process indeed surpasses top quarks as the principle irreducible background. We model the full continuum processes at NLO+PS using: We apply the same generator-level cuts as given in Eq. 3.8. We note that as events are generated using the HeavyN model file, it is necessary to also remove their diagrammatic contribution. In practice, this is done with the flag "/ nX', with nX ∈ {N 1 , N 2 , N 3 }. Alternatively, one can leave diagrams in place and either reduce active-sterile mixing elements to an incredibly small number, e.g., 10 −10 , or heavy neutrino masses to an incredibly large value, e.g., 10 10 GeV, to decouple the fields.

Electroweak (Multi-Boson) Production
Beyond the EW diboson continuum, significant but subleading backgrounds are the multiboson processes, i.e, the production of three or more EW bosons. As representative backgrounds, we consider specifically the fully resonant and mixed-resonant channels, We decay W bosons to leptons using the NWA / MadSpin in the same manner done for top quarks. Higher vector boson multiplicities, e.g., W W ZZ, suffer from phase space suppression and are ignored. Similarly, alternative weak boson permutations, e.g., W ZZ, suffer from branching fraction and coupling suppression, and hence are subleading to the processes above. For the W W ℓℓ process, we apply the generator-level cuts of Eq. (3.8).

Fake Leptons
As discussed in Sec. 1, a substantial limitation to current multilepton searches for heavy neutrinos are misidentified and mis-reconstructed/mis-measured objects [34,74], and are often neglected in theoretical studies. Though unlikely, about 1 in 10 4 QCD jets with relatively low p T are misidentified as electrons and positrons [172,173]. Similarly, it is possible for QCD jets to be misidentified as hadronically decaying τ leptons, with rates spanning about 1 in 10 2 − 10 4 as a function of jet p T [174]. Hence, assuming that a jet is mis-tagged as a charged lepton, SM processes with genuine hard lepton pairs, such as can mimic the pp → 3ℓX signal. We emphasize that not only can jets from the Born process contribute to the "fake" lepton background but also softer, wide-angle, initial-state (ISR) and final-state radiation (FSR) that are pulled into the matrix element (event file) through O(α S ) fixed order corrections and resummed corrections in the parton shower. Hence, we believe modeling the progenitors of fake leptons, where possible, to at least NLO+PS is critically necessary. To model the above processes, we use the syntax To regulate the W W j process, we impose the generator-level cuts, p j, gen.
T > 30 GeV and |η j gen. | < 3. Top quarks and W bosons are treated in the NWA and decayed to leptonic final states as above. We defer details of how fake leptons are modeled at the analysis level to Sec. 6.

Standard Model Inputs
Assuming n f = 4 massless quarks and a diagonal Cabbibo-Kobayashi-Maskawa (CKM) matrix with unit entries, we take SM inputs from the 2014 Particle Data Group [175]: We generate our model files m b , m τ = 0 in order to consistently account for hadronic decays of B hadrons and τ leptons. We nevertheless employ PDFs with variable n f . This technical inconsistency allows us to have a more realistic PDF normalization for multi-TeV m N across different colliders. PDF and α s (µ r ) evolution are extracted using the LHAPDF v6.1.6 libraries [176]. Throughout this study, we use several PDFs to best match the formal accuracy of our various collider calculations: For LO and NLO in QCD calculations with and without PS-matching, we use the NNPDF 3.1 NLO+LUXqed PDF set [177] (lhaid=324900). The PDF set features an inelastic photon PDF, i.e., photons as a parton in a hadron, matched to an elastic photon PDF, i.e., photons from hadrons, at low momentum transfer (Q 2 ) using the LUXqed formalism [157,158]. For resummed jet veto calculations at NLO+N 2 LL(veto), we use the NNPDF 3.1 N 2 LO+LUXqed PDF set (325100) to avoid double counting of O(α 2 s ) emissions. For resummed threshold calculations at N 3 LL(threshold), we use the NNPDF N 2 LO+N 2 LL(threshold) PDF set [178]. For discussions of the various PDFs and their impact on BSM collider phenomenology, including neutrino mass models, see Refs. [59,[179][180][181][182] and references therein. Similarly, while the PY8 tune we use was constructed using the NNPDF 2.3 LO PDF set, we expect the error from this mismatch to be sub-leading with respect other sources of uncertainty.
For signal and background processes, we follow the recommendation of the heavy neutrino study of Ref. [62] and set the factorization and renormalization scales to (half) the sum of all visible final-state particles' transverse energy (dynamical_scale_choice=3 in mg5amc), The impact of this scale choice on TeV-scale leptons produced via Drell-Yan is that the relative size of O(α s ) corrections to the total rate and differential cross sections are largely constant as a function of mass and the value of inclusive, leptonic observables. For further discussions, see Refs. [62,80,100] and classical references therein. Where appropriate, we estimate the impact of missing perturbative corrections to our various calculations, i.e., scale uncertainties, by varying µ f , µ r between 0.5 < µ r,f /µ 0 < 2. For PDF uncertainties, we report standard deviations across replica PDFs as prescribed by Ref. [176]. The exception to this procedure are for color-singlet processes initiated by gluon fusion, where we instead follow Ref. [66] and references therein due to the threshold resummation applied. We do not account for the perturbative uncertainties associated with the PS shower matching scale µ s . While likely important, we note that they can be brought well into control with the inclusion of additional matrix element matching [183,184].

Benchmark Heavy Neutrino Inputs
Thought this analysis, we restrict ourselves to searches for the lightest, heavy neutrino mass eigenstate, N m ′ =4 , which we denote henceforth simply as N . For discovery purposes, we also take active-sterile neutrino mixing V ℓN to be independent of light and heavy neutrino masses. This action requires care to ensure that N of arbitrarily high mass are consistently modeled in MC event generators. In accordance with the Goldstone Equivalence Theorem, for fixed active-sterile mixing in Type I-like Seesaws models, heavy neutrinos with masses much larger than the EW scale possess total decay widths that scale as Γ Tot.
Ultimately, this is due to implicitly increasing the magnitude of the Dirac neutrino Yukawa in V ℓN (to ensure small light neutrino masses), leading to N being strongly coupled, like the top quark, to the SM Higgs and the EW Goldstone bosons. For a naïve benchmark mixing of |V ℓN | ∼ 1, this cubic mass dependence is very steep: for m N 1 TeV, one obtains Γ Tot.
N /m N ∼ 1, indicating a breakdown of the particle description of N .
Even for more modest benchmark values of m N or |V ℓN |, care is still needed. For |V ℓN | ∼ 0.1, one still finds that Γ Tot. N /m N 1%, a total width of Γ Tot.
N ∼ 10 GeV means that the virtuality of a propagating, nearly on-shell N can shift naturally by δ p 2 N ∼ ±20 − 30 GeV on an event-by-event basis. This is known to broaden the characteristic kinematics of a heavy neutrino's decay products in MC simulations at √ s = 14 TeV [59] to a level that is well beyond experimental resolution, resulting in an expected experimental sensitivity to EW-and TeV-scale N that deviates greatly from true sensitivity.
On the other hand, in practice one usually obeys the relation, Γ Tot. N /m N ≪ 1, for realistic values of |V ℓN | [124]. Hence, to avoid these "large width" artifacts in our MC simulations, we develop our selection cuts and estimate our signal-significance assuming a small generator-level mixing. Results are reinterpreted for smaller mixing accordingly. For charged lepton flavor conserving scenarios, we use as our generator-level mixing, for ℓ 1 ∈ {e, µ , τ }. For charged lepton flavor-violating event samples, we use, In total, all six benchmark permutations of ℓ 1 , ℓ 2 , ℓ 3 ∈ {e, µ, τ } are considered. For more details on flavor-mixing hypotheses, see Sec. 6.2.

Heavy Neutrinos at Hadron Colliders
While production-and decay-level kinematics of heavy neutrinos are well-reported throughout the literature for collider energies spanning √ s = 2 − 100 TeV, rarely (if at all) has it been done for the express purpose of building a tailored collider analysis that remains robust against varying collider energies, let alone in the context of jet vetoes. This is a highly nontrivial task and requires one to be simultaneously aware and ignorant (inclusive) with respect to the growing transverse, hadronic activity at increasing √ s.
Therefore, in this section, we present total and differential cross sections for the production of a single heavy neutrino N and its subsequent decay to the trilepton final state, as given in Eq. 1.3 and shown in Fig. 2, across the collider energies √ s = 14, 27, and 100 TeV. We begin in Sec. 4.1, where we briefly summarize the collider production formalism that enters into our state-of-the-art predictions. In a language suitable for broader audiences, we make explicit the connections between the formal accuracy of our jet modeling and the terms that appear in the Collinear Factorization Theorem [166][167][168][169][170][171], i.e., the master equation for collider physics. In Sec. 4.2, we compare the leading heavy neutrino production channels shown in Fig. 1 and reveal a qualitative picture for heavy neutrino production that is different for the LHC and its potential successors. In Sec. 4.3, we present parton-level distributions of elementary leptonic observables that form the basis for several arguments used to build our analysis. Remaining robust against increasing hadron collider centerof-mass energies requires one to have a realistic description of QCD activity, even for the signal process. This is why show in Sec. 4.4 particle-level † distributions of both elementary and more complex kinematic observables for leptons at NLO+PS. Properties related to hadronic activity (jets) are reported in Sec. 5.
We note that due to the final-state light neutrino that appears in our pp → N ℓ → 3ℓν process, events inherently possess some degree of transverse momentum imbalance. To disentangle this from the presence of light neutrinos from τ lepton decays, we limit ourselves to the e − µ mixing / no τ -mixing scenario, with V ℓN set by Eq. 3.18. We also further restrict decays of the intermediate W boson to only e-and µ-flavored leptons. Finally, we briefly stress that for the mass range of N considered (m N > m h ) and the scales involved with the production mechanisms studied, helicity suppression of lepton number violating currents [185][186][187], a consequence of the so-called Dirac-Majorana Confusion Theorem [185] is not present. Therefore, what follows here and in Sec. 5 holds for both heavy Dirac and Majorana neutrinos produced and decayed via SM currents.

Heavy Neutrino Production Formalism
For the inclusive production and decay of a heavy neutrino N in association with X (e.g., X = ℓ or ν), we calculate total (σ) and differential (dσ) scattering rates in accordance with † Throughout this analysis, "particle-level" events, also labeled as generator-level, reconstructed events, refer to MC events that have been decayed, parton showered, hadronized, and passed through jet clustering.
the Collinear Factorization Theorem [166][167][168][169][170], given by The above states that in pp collisions, hadron-level observables (dσ) can be expressed as the product of conditional probabilities, i.e., a convolution (⊗), of process-independent parton distribution functions (PDFs) f , a quasi-universal Sudakov factor ∆, and a process-specific, parton-level ik → N X hard scattering observable dσ, which occurs at the hard scattering scale Q = m N X ∼ m N , all up to suppressed corrections that scale as powers of the hadronic, or non-perturbative (NP), scale Λ NP ∼ 1 − 2 GeV.
More precisely, f i/p (ξ, µ f ) represents the likelihood of observing a massless parton i ∈ {q, q, g, γ}, where q ∈ {u, c, d, s}, carrying a transverse momentum p i T below a collinear factorization scale µ f and a longitudinal momentum p i z = ξP z away from a proton p with momentum P ≈ ( √ s/2, 0, 0, ± √ s/2). Via the DGLAP renormalization group evolution equations [188][189][190], f i/p accounts for (resums) all collinear j → i + j ′ QCD and QED emissions leading to parton i such that p i, j ′ T < µ f , i.e., f i/p accounts for all initial-state gluons, photons, and quarks with p T < µ f emitted in association with the (N X) system. Radiations with p j ′ T > µ f are included through appropriate O(α m s α n ) perturbative corrections to the hard matrix element dσ. Inclusion of QED DGLAP evolution [155,157,158,191] that is matched at µ f → m p ∼ 1 GeV to a proton's elastic photon PDF implies [155] the photon PDF f γ/p is valid for all photon virtualities. This is needed to correctly model inclusive N ℓ ± production from W γ fusion [59,62]. Together, partons i and k from PDFs f i/p and f k/p generate what is labeled as the parton scattering, or partonic, scale √ŝ = √ ξ 1 ξ 2 s ≥ Q. Intuitively, the Sudakov factor ∆(z) = δ(1−z)+O(α s )+O(α), through exponentiation and resummation [192][193][194], acts to "dress" (in the quantum field theoretic sense) bare partons with collinear and/or soft, i.e., unresolved, QCD and QED radiation, rendering them more inline with physical quantities. (In practice and reality, arbitrarily infrared emissions are regulated by hadronization.) More technically, the Sudakov factor accounts for gluons (and photons) carrying a momentum fraction (1 − z) emitted prior to the ik → N X hard process (σ ββ ′ ). Here, z = Q 2 /ŝ is a dynamic measure of how much of the total partonic energy is carried by the hard (ik) process. ∆ ββ ′ ik is quasi-process-independent in the sense that it is sensitive to the color and Lorentz structure of the incoming partons, e.g., ββ ′ = F F for (qq ′ )-or AA for (gg)-annihilation into a color singlet final state, but insensitive to the hard process itself. The precise organization and exponentiation of O(α s ) (and higher order) terms in ∆ ββ ′ ik (z) leads to the all orders resummation of logarithms associated with collinear and/or soft ISR and FSR. This includes terms associated with jet vetoes that are present in differential distributions but cancel for the inclusive cross section [92,168,170,[195][196][197][198][199][200]. To recover fixed order (FO) predictions from Eq. 4.1, one simply neglects all terms in ∆(z) beyond the δ(1 − z) kernel. Notably, collinear ISR and FSR are well-modeled by Sudakov form factors as implemented in modern parton showers (PS) since they formally resum corresponding terms to at least leading logarithmic (LL) accuracy, when expanding the color algebra in powers of 1/N c [201].
For the n-body ik → N X hard process, the parton-level observablesσ and dσ, for unpolarized hadron collisions, are obtained from the expressions, Here, s j = 1/2 are the helicity degrees of freedom (dof) for massless parton j = i, k; N j c are similarly the SU(3) c color dof; M is the Lorentz-invariant matrix element calculable via perturbative methods; and dP S n is the separately Lorentz-invariant n-body phase space. For NLO in QCD corrections, we employ the MC@NLO [152] formalism as implemented in mg5amc in order to insure real radiation common to the PDFs f , Sudakov factor ∆, and dσ beyond LO are appropriately subtracted, up to our claimed NLO precision.
Beyond perturbative corrections to matrix elements in dσ, to all orders, the impact of collinear, initial-state QCD and QED emissions on the normalization of the total, hadronic cross section (σ) are, by construction, already included in the DGLAP evolution of PDFs † . This is the reason for subtraction of O(α s ) terms from the PDF when dσ is known at NLO in QCD. By unitarity, collinear FSR does not readily impact the normalization. Qualitatively, however, both ISR and FSR can significantly alter the shape of hadron-level differential distributions (dσ) and should not be neglected. Independent of jet vetoes, neglecting the parton shower entirely and considering searches only at the (unphysical) partonic level has a substantial impact on selection cut efficiencies for Seesaw particles [34,47,49].
The impact of soft/threshold logarithms [205][206][207] contained in ∆, which originate from considering additional partonic phase space contributions [208], can sizably increase the total cross section normalization with respect to Born-level predictions for certain kinematic regimes [209] and initial-state color configurations [210]. (Though these resummed contributions are not entirely independent of FO contributions.) For new physics processes, this is still the case after removing potential double counting of soft-collinear contributions in PDFs [66,[178][179][180][181][182]. For the GF production process in Fig. 1(b), and exploiting the renormalization group (RG) evolution properties of the Factorization Theorem in momentum † For clarity, the normalization increase of the trilepton process in Fig. 2 at NLO in QCD is almost entirely driven by the virtual correction and results from a combination of QCD color factors, a loop phase space suppression factor, and a π 2 term due DY to being a space-like process [202][203][204]. The typical correction of +20 to +40% at NLO does not original from ISR; see, e.g., Ref. [80] and references therein.
space [148,150,151], the N 3 LL threshold corrections can be obtained in a straightforward manner [66]. The procedure ostensibly requires replacing the FO expression for the Sudakov factor ∆ FO ≈ δ(1 − z) with the resummed / RG-improved expression at N 3 LL(threshold), ∆ Res. , while adding and subtracting the necessary terms to coax numerical stability. For NLO+N 2 LL jet veto corrections, the procedure is more involved but nonetheless straightforward [96]. As for pure fixed order corrections, it requires subtraction of collinear splittings from the PDF f . However, instead of simply augmenting α s -subtracted PDFs with α scorrected matrix elements, one introduces corrections that imbues incoming / initial-state partons with transverse momentum. These corrections act as an intermediate term connecting (a) the PDFs evolved up to the collinear factorization µ f to (b) the p T of the leading QCD emission in the matrix element. (In some sense, the result is akin working with transverse momentum dependent (TMD) PDFs that are matched to the α s -corrected matrix element.) Power-suppressed, non-perturbative terms of the order O Λ n NP /Q n+2 , i.e., higher twister terms in the operator product expansion [166,167,169,211], beyond hadronization are largely neglected in this study. For jet vetoes, this is a potential source of sizable theoretical uncertainty. More precisely, the momenta of reconstructed jets in hadron collisions, and hence jet vetoes cross sections, receive corrections from Glauber exchanges, e.g., double parton scattering/multiple parton interactions/underlying event (UE). Such corrections are beyond the theorem's formal accuracy nor are presently known to fully factorize; for further details, see Refs. [168-171, 200, 212, 213]. While we neglect their impact, we do believe it is possible to reliably parametrize UE, particularly at larger √ s, as demonstrated in Ref. [214]. A dedicated investigation into the impact of UE on jet vetoes applied to new physics searches is deferred to future studies. For other related details, see Sec. 5.

Factorization of Heavy Neutrino Mixing from Collider Observables
For the resonant production of a single heavy neutrino in pp collisions, the active-sterile mixing element V ℓN enters the scattering matrix elements as an overall multiplicative factor. Hence, the (squared) norm of the active-sterile mixing, |V ℓN | 2 , factors out at the squared matrix element, and therefore factors out of at the hadronic level. This gives rise to the so-called "bare cross section" dσ 0 , defined, for example for pp → N ℓ, as [46], as Subsequently, semi-flavor-model independent prediction for both Dirac and Majorana neutrinos can be constructed at the N production level. |V ℓN | 2 also enters into the partial decay widths of N , and hence its total decay width. As this in turn appears in the Breit-Wigner propagator for N , the factorization of neutrino mixing elements from decay-level heavy neutrino observables is slight more subtle. The analogous "bare cross section" for resonant production and decay of N , in for example Notably, both factorizations hold identically to all orders in QCD [62,80], so long as the individual decay products of N remain color neutral; for hadronic decays of N , the factorization holds at least to NLO. Trivially, dσ 0 can be obtained numerically by setting |V ℓN | = 1 in dσ or by a rescaling for nonzero |V ℓN |.
For multiple interfering heavy neutrino mass eigenstates, the above relations cannot be generically applied. This is particularly the case for pseudo-Dirac neutrinos, where relative CP phases in the V ℓN lead to significant destructive interference for L-violating processes. However, as discussed in Sec. 2.1, the mass splitting for pseudo-Dirac neutrino pairs in a realistic ISS should be vanishingly small on the collider scale. Hence, pseudo-Dirac pairs can be justifiably modeled as a single Dirac neutrino at the collider scale.

Heavy Neutrino Production at the LHC and Beyond
To begin quantifying the discovery potential of trilepton searches for heavy neutrinos at pp colliders, we compare the leading cross sections for a single resonant N produced through EW bosons. Using Eq. 4.1 and the methodology outlined in Sec. 3, we plot in Fig. 3 the bare production cross section σ 0 (as defined in Eq. 4.5), for N produced though the CC DY and NC DY channels at NLO in QCD, as shown in Fig. 1(a) and given by, at the Born level. These processes, especially the CC DY, are to serve as baselines for subsequent discussions. We also plot the N 3 LL GF mechanism as shown in Fig. 1 as well as the W γ VBF channel at NLO in QCD, shown in Fig. 1(c), and given by where σ N j LO+N k LL is the cross section evaluated at N j LO+N k LL in QCD and σ LO is the LO/Born cross section. For the 14 TeV LHC, the nontrivial interplay between matrix elements and PDFs for the above processes has been reported in detail elsewhere [37,38,41,58,59,62,66] and will not be discussed further. We note simply that for the mass range m N = 150 − 1000 GeV, the bare cross sections and their relative uncertainties span roughly  A detailed summary of heavy N production rates, their residual scale dependencies, and their QCD K-factors, for representative m N and √ s is given in Tab. 1. In all cases, the largest uncertainties are associated with the lowest mass points and quickly moderate to smaller levels. This falloff is particularly notable for the DY channels. The sensitivity at the lowest (m N / √ s) is attributed to the opening of qg-and γg-initiated processes, and hence sensitivity to the gluon PDF, which slopes most significantly at low Bjorken-x. Moreover, at LO, the W γ fusion process only possesses µ f dependence (through the PDF); µ r dependence arises through α s first at NLO. The GF channel, despite its high accuracy, possesses a large factorization scale uncertainty due to missing real radiation terms that would otherwise stabilize the prediction [66,149]. Qualitatively, one sees a number of important changes in the relative importance of heavy neutrino production mechanisms when going from √ s = 14 TeV to 27, 50, or 100 TeV. While in all cases, the VBF process overtakes the CC DY process at a largely static m N ∼ 900 GeV, the dominance of GF varies due to the surge in the gg luminosity for diminishing (m N / √ s). The constant cross over for the CC DY and VBF channel is likely tied to the fact that both are driven by valence quark-sea parton scattering, and therefore receive the same growth in parton luminosity as collider energy increases. At 14 TeV, GF is sub-leading for all m N under consideration but comparable to the NC DY process. Technically, the NC DY and GF should be summed coherently as the GF channel is formally a separately gauge invariant O(α 2 s ) correction to the NC DY mode [41]. (Similarly, the W γ channel is a separately gauge invariant EW correction to the inclusive CC DY process [59].) Doing so renders it much closer to the (CC DY+VBF) process at 14 TeV [66]. At 27 (50) TeV, GF emerges as the leading channel for m N ≈ 450−900 (300−1200) GeV, beyond which VBF is dominant. Process pp → N νX process, few (if any) dedicated collider studies exits for m N > m h . However, this regretful absence is partly due to the process' larger SM backgrounds.
The importance of the VBF channel cannot be overemphasized: across the four colliders, for m N = 1 (2) [3] {10} TeV, the bare cross sections reach about σ 0 ∼ 1 − 6 fb. In light of the proposed O(3 − 30) ab −1 datasets being earnestly discussed by the community [29][30][31]33], one expects that active-sterile neutrino mixing on the order of |V ℓN | 2 ∼ 10 −4 − 10 −2 for this mass range can be probed by the VBF channel alone. In the context of Type I-like Seesaw models, these rates also suggest the first direct sensitivity to O(10) TeV heavy neutrino masses as colliders, and compels us to focus on the discovery prospect of the inclusive (CC DY+VBF) pp → N ℓX signature for the remainder of this study.

Parton-Level Kinematics and Observables at LO
In this section, we present at √ s = 14, 27, and 100 TeV, the parton-level kinematics of the trilepton process for a single heavy N via the CC DY mechanism, at LO in QCD. The purpose of this section is to establish baseline kinematical features inherent to the qq ′ → N ℓ → 3ℓν hard process and explore the dependence, or lack thereof, on increasing collider energies. Observations established here ultimately form the basis for several conclusions we make at the particle-level in Sec. 4.4 as well as justifications for analysis selection criteria.
As we are working momentarily at the partonic level, we have access to the standard HEP MC particle identification codes (PID) of final-state particles and their resonant parent particle in the outputted mg5amc LHE files. Therefore, for each charged lepton in the following figures, the subscript label X ∈ {N, W, ν} denotes the particle with which the charged lepton was produced: ℓ N (solid) denotes the primary charged lepton produced in association with N ; ℓ W (dash) represents the charged lepton from the N decay and produced with an on-shell W boson; and ℓ ν (dot) is the charged lepton from the leptonic W decay. Notably, the recording of the PID and 4-momentum is true even for the intermediate N , which is not modeled in this section with the NWA (see Sec. 3.2). In MC event generation with mg5amc, intermediate particles that go on-shell are recorded in the event file independent of the NWA being applied. As the widths of heavy neutrinos used in MC generation are engineered to be very small due to small active-sterile mixing (see Sec. 3.5), in accordance with realistic neutrino mass models. We find nearly all ( 98%) neutrino events generated at NLO+PS contain an on-shell N with its momentum recorded to file.
Briefly, we note that we omit kinematics associated with the W γ fusion production as they have been studied systematically for various √ s elsewhere [59]. In addition, the benchmark (pseudo)rapidities † for nominal tracker coverage of charged leptons at future collider studies are set to y, η = ±2.5 [27]. Hence, throughout this section, guide lines at these values are inserted in all (pseudo)rapidity plots.
We begin with a look at the kinematics associated with the hard (N ℓ)-system itself. In Fig. 4 we plot the normalized invariant mass distribution of the (N ℓ)-system for (a) m N = 150 and (b) 450 GeV, at √ s = 14 (solid), 27 (dash), and 100 (dot) TeV. Immediately, one observes that for a fixed m N the shape of the invariant curves are largely the same despite the large differences in collider energies; only minor broadening is observed for increasing √ s. The reason is that up to kinematic factors related to phase space and angular momentum (e.g., s-vs p-wave), the dominant contribution of any matrix element is in the neighborhood of its poles, i.e., when particles are resonant and go onto their mass shell, and is (mostly) independent of collider configuration. This is seen explicitly in the † For massive 4-momentum p µ , rapidity is defined as In the massless limit, y simplifies to pseudorapidity, y where θ is the polar angle with respect to the beam/ẑ-axis.
with Q m N denoting the scale of the hard process, which in this case is equal to the invariant mass of the (N ℓ)-system. In the Q, m N ≫ M W limit, one sees clearly the kinematic threshold at Q = m N , below which the process is kinematically forbidden (by momentum conservation). The kinematical factor of (2Q 2 + m 2 N ) acts as a brief "turnon" function before the ∼ 1/Q 8 suppression sinks the rate for larger invariant masses. Numerically, the maximum occurs as Q ∼ (4/3) × m N , and indicates that the hard process scale Q is of the same order as the threshold scale m N , independent of collider energy. Now, as the 2 → 2 process pp → N ℓ N occurs at a characteristic invariant mass of Q ∼ m N , the 4-momenta of N and ℓ N should also carry characteristic values of this order.
To demonstrate this, we show in Fig. 5, the normalized (a,b) transverse momentum (p N T ) and (c,d) rapidity (y N ) distributions of N for (a,c) m N = 150 and (b,d) 450 GeV, at √ s = 14, 27, and 100 TeV. Like for the invariant mass of the (N ℓ)-system, for a fixed m N , one observes that the shape of the heavy neutrino p T distribution, which scales as is largely constant despite the large differences in collider energies, confirming the conjecture at LO. On the other hand, the y N distribution broadens significantly with increasing collider energy. This latter feature can be attributed to the fact that, at LO, increasing √ s translates directly to increasing the range of longitudinal momentum that initial-state partons may hold, and hence the range of longitudinal boosts imparted on the (N ℓ N )-system and its descendants. This suggests that observables that depend on the longitudinal momentum of N or the (N ℓ N )-system will tend to broaden with increasing √ s. As a function of increasing    Figure 5. distributions of the charged lepton produced in association with the heavy neutrino (ℓ N ) mirror those of N and need not be discussed further. The η distributions for ℓ N are slightly narrower than the y N curve due to ℓ N 's (approximately) massless nature.
Keeping to Fig. 6, for both low and high mass N , one sees that the distribution shapes are largely independent of collider energy. This follows from the fact that the p T of ℓ W and ℓ N have characteristic values due to the intermediate resonant structure of the N → ℓ W W → ℓ W ℓ ν ν decay. Consequently, up to relative transverse boosts of N , the transverse momenta of ℓ W and ℓ N scale as, where in the last approximation we take the (M W /m N ) 2 ≪ 1 limit. As we have argued, however, for the inclusive production of N , p N T is largely stable against varying beam energy and therefore we expect and observe this independence to be inherited.
For the η distributions in Fig. 7, a number of notable features are observed: The first is that for a fixed √ s and m N , the three pseudorapidity curves overlap considerably. This follows from the combined circumstance that each charged lepton possesses a comparable transverse momentum, one that scales like m N , and that all three charged leptons obtain    their relative longitudinal momentum from the same source, namely the hard (N ℓ N )-system. As such, the distributions are indistinguishable and cannot be readily used to help identify a particular lepton's origin. Second, for the narrowing of the η distribution with increasing m N , we reiterate that this simply reflects the relative increase of p T observed in Fig. 6.
Regarding the √ s dependence, again one observes a broadening with increasing collider energy, which is inherited from the parent (N ℓ N )-system. That said, what is important to stress is the degree of the migration of events outside the |η| = 2.5 boundary. For m N = 150 GeV, we see that the fraction (f ) for each lepton to fall outside this detector boundary grows from about f 20%, to 30%, to 40% as one goes from √ s = 14 TeV to 27 TeV and 100 TeV. Neglecting correlation, this means that only about ε = (1 − f ) 3 50% (35%) [25%] of trilepton heavy neutrino events fall within the benchmark fiducial coverage at √ s = 14 (27) [100] TeV. For m N = 450 GeV, the situation is marginally better with ε 85% (60%) [35%], but still discouraging. For even larger m N , one can be more hopeful, but at the cost of a smaller cross section (for a fixed active-sterile mixing).
All-in-all, as evident in Fig. 7, extending tracker coverage to at least |η| = 3 − 4 for future collider detector experiments would be vital, if not necessary, to ensure a healthy search program for anomalous production of charged leptons. When these distributions are taken together, a fuller picture starts to emerge at LO: for symmetric, anti-collinear beams as we consider, as one increases √ s, the incoming initialstate partons entering the qq ′ → N ℓ hard process are allowed to carry larger longitudinal momenta along the beam axis. As a result, the (N ℓ N )-system itself is allowed to carry a larger longitudinal boost, and hence the observed broadening of the y N spectra. On the other hand, to an excellent approximation, the net p T if the initial proton-proton system is zero and does not change as a function of collider energy. This suggest that the p T of the (N ℓ N )-system does not depend on √ s, which is essentially true at LO, and is reflected in the p T distributions of N , its descendants, and ℓ N . At the leading log / parton shower level, however, this argument of course begins to break down: there is soft/collinear ISR that imbues the (N ℓ N )-system with nonzero p T that does indeed depend, to some extent, on the phase space made available by the collider center of mass energy. This radiation, however, dominantly contributes to the region of phase space where p That is to say, the p T of the (N ℓ N )-system is typically small compared to the invariant mass of the (N ℓ N )-system [166,167]. And as demonstrated in Fig. 4, the invariant mass of the (N ℓ N )-system is driven by the matrix element, not by √ s. So while we do expect the p T of N , its descendant leptons, and ℓ N to change with increasing √ s, this change will be small compared to the hard scale Q ∼ m N , which sets the initial p T scale of N and ℓ N .
Interestingly, this then suggests that kinematic observables derived strictly from the p T of ℓ N , ℓ W , and ℓ ν will inherit this robustness against changing √ s. We now turn to particle-level kinematics at NLO+PS(LL) to see how far this holds.

Particle-Level Kinematics and Observables at NLO+PS(LL)
In this section, we investigate particle-level observables for the same trilepton process in Eq. 4.27, assuming the same representative m N and √ s, but instead at NLO+PS. We do this by first highlighting the key distinctions between parton-and particle-level objects (these are passed through parton showers, hadronization, and jet clustering algorithms). We then point to a subtlety in MC generation at NLO+PS that enables us to establish a definitive link between parton-and particle-level kinematics. Building on this correspondence, we are able to see, using this significantly more realistic modeling, how robust the kinematics of N and its descendants from the process in Eq. 4.27 are against variable collider energy. We then demonstrate the existence of an entire class of particle-level observables that exhibits only a weak dependence on √ s. This is a main result of this study.
For several reasons operating at the parton shower-level with hadronization modeling is crucial, particularly for this study: While parton-level kinematics provide important insight and intuition to underlying processes, not all kinematic observables built at the parton level are physical in hadronic environments [166,167]. Moreover, detector experiments do not observe partons or bare charged leptons: at macroscopic distances, high-p T partons and charged leptons are dressed with collinear and soft QCD and/or QED radiations. This has a fundamental impact on both exclusive observables, such as jet and lepton multiplicity, as well as inclusive ones, such as the transverse momentum of the (N ℓ N )-system. In addition, when the hard scattering process and beam remnant separate beyond a distance of d NP ∼ 1/Λ NP ∼ 1 fm, color-connection and hadronization must be taken into account. This gives rise to a cornucopia of (relatively) low-p T hadrons and such contributions play an important role in this work: Depending on the largeness of the jet radius parameter R, "out-of-cone" and "into-cone" emission of hadrons can shift jet p T by up to O(10) GeV [212]. Likewise, light neutrinos originating from weak decays of hadrons can cumulatively shift the direction and magnitude of the missing transverse momentum vector by O(Λ NP /decay). Not to mention, hadronic decays of τ leptons (τ h ) require modeling and experimental tuning as their branching rates cannot be derived from first principles.
Despite the activity of a hadron collider environment, it is possible to make a more or less concrete mapping to parton-level objects, though not necessarily their precise identities. As done in detector experiments, one must cluster (sum) the momenta of all like-objects into composite objects (jets) according to some procedure (a sequential jet clustering algorithm) that is valid at all orders of perturbation theory (is infrared and collinear-safe). With this additional layer of abstraction, a many-body environment is simplified to a few-body system with correspondence to the partonic event. This procedure, though, is not entirely free of ambiguities. A consequence of working at this reconstructed level is the obfuscation of amplitude-/generator-level information about a particle's particular origin. For example: in Figs. 6 and 7, the degree of overlap and similarity of the p T and η curves for ℓ N , ℓ W , and ℓ ν is worsened due to recoils against electromagnetic FSR and QCD ISR. In effect, the exact lineage of the charged leptons are anonymized at the reconstruction level.
A quirk exist, however, in the present case that allows us to establish a stronger link between parton-level and particle-level heavy neutrino events. As described in Sec. 4.3, for event generation with mg5amc, intermediate particles that go on-shell are recorded in the event file independent of the NWA. This is true even at NLO. Moreover, the 4-momentum and PID of intermediate, on-shell particles from the original qq ′ → N ℓ → ℓℓW hard process are preserved throughout the MG5aMC@NLO+PY8+MA5 chain, thereby allowing us to build kinematic distributions of N at this accuracy.
Subsequently, we show in Fig. 8  one sees very little changes in any of the curves and, crucially, that the insensitivity to √ s is retained. This is due to the fact that QCD corrections for DY-like processes are dominated by positive virtual corrections [202][203][204]. This further indicates that the p N T and y N distributions (and the analogous ones for ℓ N ) exhibit Born-like features at NLO+PS  Figure 8. Same as Fig. 5 but at NLO+PS accuracy with particle-level reconstruction.
but with a possible enhancement at large p T owing to the opening of partonic channels, e.g., Compton scattering gq → N ℓq ′ with p q ′ T > µ f . However, following Ref. [62], our choice of factorization and renormalization scales, as given in Eq. 3.16, are intentionally selected to minimize the impact of such channels on distribution shapes. Thus, such enhancements will be minor, and one may infer that observables directly related to the kinematics of N , e.g., p T of its decay products, remain mostly unchanged relative to LO predictions.
Focusing on the decay products of N , we plot in Figs Due to QED parton shower effects and leptonic decays of hadrons, more than three charged leptons may be present in each event. Hence, here and for the remainder of the text, leptons ℓ k with k = 1, . . . , are ordered such that p ℓ k T > p ℓ k+1 T . Despite the potential increase in charged lepton multiplicity, one sees in Fig. 9 characteristic structures emerge in the p T curves. This indicates that the three leading charged leptons largely originate from heavy, intermediate resonances, and are not low-energy, continuum contributions. At m N = 150 GeV, ℓ 1,...,3 cannot be reliably or uniquely associated with ℓ N , ℓ W , or ℓ ν , which is expected as the LO curves in Fig. 6 indicate that the p T are all highly comparable for fixed m N . For heavier N , the association is slightly stronger, with one observing that p ℓ 1 T ∼ m N /2 and p ℓ 2 T ∼ m N /4. This suggests that the leading charged leptons may be more readily identified with ℓ W and ℓ ν , respectively, but still with large uncertainty. In the pseudorapidity curves of Fig. 10, no discerning feature is present to discriminate or associate ℓ 1,...,3 with ℓ N , ℓ W , or ℓ ν . Even the slightly narrower/taller distribution for η ℓ N seen at LO is washed out at  Figure 9. Same as Fig. 6 but at NLO+PS accuracy with particle-level reconstruction.
NLO+PS. In short, upon discovery of an anomalous trilepton signature consistent with heavy neutrino production and decay, it will be difficult to readily identify the event's three leading charged leptons with ℓ N , ℓ W , or ℓ ν . To do so, one must resort to more sophisticated techniques, such as the Matrix Element Method, more complex observables, or additional hypotheses, such as lepton number/flavor conservation or violation.
Regardless of the ability to readily associate ℓ 1,...,3 to ℓ N,W,ν , with respect to the dependence on collider energy, both classes of the p T and y/η distributions reflect the behavior  Figure 10. Same as Fig. 7 but at NLO+PS accuracy with particle-level reconstruction.
observed in the LO distributions. Subsequently, for a fixed heavy neutrino mass, one can conclude that since the shapes of the p T distribution for N remain effectively insensitive to changes of √ s, so too do the distributions for ℓ 1 , ℓ 2 , and ℓ 3 at NLO+PS. We now investigate how far this robustness against varying collider energy holds and consider more complex observables that are derived from the p T of the leading charged leptons.
In addition to the three charged leptons, the collider process of Eq. 4.27 also contains final-state light neutrinos originating from the N → ℓW → 2ℓν decay. As mentioned above,  Figure 11.
where the summation runs over all visible final states regardless of their hardness (p T ). Disregarding objects below a particular p T threshold can introduce distortions in the MET distribution of that same order and worsen the perturbative stability [100]. In Fig. 11, we plot the particle-level MET for (a) m N = 150 and (b) 450 GeV, at √ s = 14, 27, and 100 TeV. Immediately, one sees that MET, a quantity built from both leptonic and hadronic transverse momenta, more or less aligns with naïve 1 → 2 → 3 kinematics for N decays, which stipulate that the light neutrino p T and E scale as (115)  Like the p ℓ X T distributions in Fig. 9, the structured behavior indicates that the distributions are driven by heavy, intermediate, resonant states. In contrast to the p T of the leading charged leptons, we see a stronger collider energy dependence. The MET peak drops 10to-20% from its maximum as one goes from 14 TeV to 100 TeV, for both m N = 150 and 450 GeV. The location of the maximum, however, does not noticeably change. This broadening can be attributed to two effects: (i) At higher collider energies, more low-p T hadrons are generated due to the larger available phase space. In turn, more hadrons can undergo leptonic weak decays, which in turns generates more light neutrinos. (ii) Deviations from Eq. 4.34 also originate from neutrinos carrying momentum in the longitudinal direction, which worsens with increasing √ s due to increasingly forward valence quark-sea antiquark  annihilations. This can be seen in the rightward shifts of the MET tail. So while the p T of the light neutrino appearing in the N → ℓW → 2ℓν decay may remain robust against √ s, its proxy (MET) is slightly less robust due to hadronic contamination.
One should caution on the strength of the conclusion drawn from Fig. 11. Here we only show particle-level MET and do not (yet) account for detector effects. In reality, higher collider energies give rise to a higher multiplicity of jets, which can cause additional MET through momentum mis-measurement and finite fiducial coverage.
To investigate the collider energy dependence on global leptonic activity, we consider two representative measures: the (exclusive) scalar sum of leptonic p T , S T , defined as with the sum running only over the three leading charged leptons in the event, as well as the invariant mass of the same three charged leptons: (4.36) In Fig. 12, we show S T for (a) m N = 150 GeV and (b) 450 GeV, and likewise m 3ℓ , respectively, in (c) and (d). As anticipated, only a slight broadening of S T and m 3ℓ occurs with increasing collider energies. We observe that all four peaks only reduce 5 to 10% from their maximum as one goes to 27 or 100 TeV. Like MET, we see that the location of the peak remains unchanged. The weaker dependence on √ s is attributed to the neglect of charged leptons with p T < p ℓ 3 T . As the hard pp → N ℓ → 3ℓX process only contains three charged leptons in the final state, additional charged leptons in the event must necessarily originate from the decays of the low-p T hadrons in the beam remnant and possibly electromagnetic ISR/FSR. As we saw in Fig. 11, these additional contributions have a collider energy dependence, and therefore neglecting their contribution to S T and m 3ℓ helps fortify these observables' insensitivity to √ s. It is noteworthy that both observable peaks near S T , m 3ℓ ∼ m N . This is a somewhat a fortuitous accident and results from the fact that the momenta of the three charged leptons, as observed in Fig. 6, scale as For a heavy neutrino with masses that we consider, this means for S T we have Whereas S T slightly undershoots m N , we see that m 3ℓ lurches above it due to the additional longitudinal momentum that goes into m 3ℓ , which also explains why the m 3ℓ distributions are broader than the S T curves. Otherwise, there is little difference between the observables. However, due to this broadening, we prefer to build and limit our signal analysis to quantities constructed only from transverse momenta. We stress though that while S T and m 3ℓ are numerically comparable to m N for the CC DY process, this does not hold for all production mechanisms. The W γ fusion channel, for example, kinematically favors lower p ℓ N T due to an initial γ → ℓℓ * splitting [59]. Nonetheless, taking for example the extreme limit of (p ℓ N T /m N ) → 0, one still observes that S VBF T 3m N /4. Hence, the more broadly applicable statement is that S T scales proportionally with m N in a somewhat universal fashion, a fact that we will exploit in our signal analysis.
From both a search-strategy and a properties-measurement perspective, being able to reconstruct the mass of N is highly desirable. The capability, however, is complicated by the light neutrino in the N → 2ℓν decay chain, which makes it impossible a priori to reconstruct reconstruct N 's 4-momentum, and hence its invariant mass, from the charged leptons and MET momentum vectors alone in pp collisions. And as seen in Fig. 12, leptonic observables such as S T and m 3ℓ do not appear sufficiently robust estimates of heavy neutrino masses due to their broadness. Though interestingly, like the top quark and Higgs boson, it is possible to build a more reliable measure of the invariant mass of N . We do this by exploiting simultaneously the 1 → 2 → 3 decay structure of N and that the intermediate W boson is largely on its mass shell for the range of m N we consideration. For 1 → 3 and 1 → 4-body decays, a class of transverse mass observes exist [215][216][217] that are essentially multi-body extensions of the canonical, two-body transverse mass variable used in W boson mass measurements [215]. Differences in the observables are based on the multiplicity of final-state neutrinos and whether decays are sequential or in parallel.  For our N → ℓ W W → ℓ W ℓ ν ν configuration, two † possible mass-estimates can serve our † Using comparable assumptions, techniques pioneered for leptonic and dileptonic decays of top quark pairs [218][219][220][221] may also be applicable to heavy neutrinos. Such investigations are encouraged.
purpose: the multi-body transverse mass m M T (labeled as M ′ T,W W in Ref. [217]), and the multi-body cluster mass m Cl (labeled as M C,W W in Ref. [217]), Considering the assumptions that go into constructing these quantities, it not at all obvious whether one observable is (or should be) a better estimate of m N . Therefore, we compute both: At the particle-level, charged leptons ranked according to p T cannot readily be associated with ℓ N , ℓ W , or ℓ ν (see Fig. 9). By electric charge conservation, however, we know that the sum of the three's electric charges must add to Q 3ℓ = 3 ℓ Q ℓ = ±1 and that the charges of ℓ W and ℓ ν must cancel. Hence, without guidance, there are in principle two (four) permutations of ℓ 1,...,3 that one could use to build either mass variable for Dirac (Majorana) N . One choice should concentrate near the true invariant mass of N whereas the others should follow a continuum distribution. Due to the low multiplicity of choices, we (in an unsophisticated manner) brute force a guess-and-test determination. Specifically, for a particular mass hypothesis (m Hypo. To help gauge whether one transverse mass better approximates m N , we take advantage of our accessibility to generator-level information on N and set m Hypo. N , a qualitative difference is noticeable. Whereasm M T appears to better describe high mass N , i.e., give a sharper, narrower distribution,m Cl is found to do better for lower N masses. We do not investigate this difference further, but note that such behavior may be useful input in constructing multivariate analyses. In light of the better performance for higher mass N as well as our goal of quantifying the discovery potential of EW-and TeV-scale N at higher energy pp colliders, we focus on employing the multibody transverse mass m M T for the remainder of this study.

Summary
In this section we have reported in exhaustive detail production-and decay-level properties of high-mass (m N > m h ) heavy neutrinos produced in hadron colliders, at various perturbative accuracies and collider energy of mass energies. Due to the amount of detail, it is appropriate to briefly summarize the immediate findings: We report that the leading production mechanisms of heavy neutrinos across a range of masses varies enormously as a function of √ s. While the W γ fusion process inevitably dominates for the heaviest mass scales, the largely studied CC DY process becomes subleading to neutral current processes at collider energies immediately beyond those reached at the LHC. For a fixed heavy neutrino mass m N , we also report the existence of an entire class of observables, namely exclusive quantities built (dominantly) from the transverse momenta of the leading charged leptons, whose distribution shapes display only a weak sensitivity to changes in collider energy.
We now turn to the jet activity in heavy neutrino production at pp colliders.

Heavy Neutrinos and Dynamic Jet Vetoes
We now investigate the hadronic activity in heavy neutrino events when they are produced through the CC DY and W γ mechanisms. We consider this for a range of neutrino masses and hadron collider energies, as well as for representative SM backgrounds to the inclusive collider signature. It is important to reiterate that here and throughout the text we mean "inclusive" with respect to hadronic activity and, potentially, additional charged leptons that may appear in an event. The goal of this section is to quantify the amount of hadronic activity in these various processes and to establish to what extent the relative amount of hadronic and leptonic activities in these processes changes with √ s. As noted in Sec. 4.4, for a fixed heavy neutrino masses, the amount of leptonic activity, as measured by charged lepton p T or exclusive S T , does not change appreciably with √ s. The main result of this section is to show that discriminating events according to the relative amounts of leptonic and hadronic activities, on an event-by-event basis, is a powerful means to improve the signal-to-background ratio (S/B) in multilepton searches for heavy N . We reach this conclusion in the following manner: In Sec. 5.1, we discuss the characteristic behavior of QCD activity in the CC DY and W γ signal processes as a function of m N and √ s. We present how jet vetoes, which exploit this information, can be extended to account for the relative amount of leptonic activity on an event-by-event basis, i.e., a dynamic jet veto. In Sec. 5.2, we look into how a dynamic jet veto impacts leading backgrounds. We briefly discuss alternative, complementary measures of hadronic and leptonic activities that may be fruitful for future studies in Sec. 5.3. Lastly, in Sec. 5.4, we examine the uncertainties associated with static and dynamic jet vetoes at the parton shower level.

Signal Processes and Jet Vetoes at NLO+NNLL
As shown in Fig. 3, the inclusive trilepton process of Eq. 5.1 is driven by the CC DY mechanism, i.e., quark-antiquark annihilation, at lower neutrino masses and by W γ fusion for larger masses, These are both color-singlet processes where initial-state partons scatter into a colorless final state. This often leads to the naïve (and incorrect) suggestion that there is no QCD radiation in hadronic production of heavy neutrinos, especially at leading/lowest order. Protons, of course, are color-singlet states. Hence, for DY-like process, whatever remains after removing the initial-state quark or antiquark (a 3 or 3 under QCD) must also be collectively charged (as a 3 or 3, respectively) under QCD due to color conservation. As a result, initial partons and the beam remnant are connected via color fields that spontaneously give rise to jets that are characteristically collinear with the beam axis. For the W γ process, the argument is less subtle: initial-state W bosons are generated perturbatively in pp collisions from q → W q ′ splitting. As a result, the W γ → N ℓ process is typically associated with a forward jet possessing p j T ∼ M W /2 that remains spatially close and color-connected to the rest of its parent proton. (A subleading, high-|η| jet with p T ≪ M W /2 in backward direction is also usually present.) Beyond LO, this behavior remains much the same for both SM processes [163,[202][203][204][222][223][224][225] and heavy neutrino production [61,62,80].
Regardless of production mode, color conservation and mass scales lead to the fact that heavy neutrinos are accompanied by jets that are predominantly forward (large |η|) and/or soft (low p T ). This is in contrast to QCD processes, such as top quark production, where jets are typically central (small |η|) and hard (high p T ). This qualitative difference provides the rationale for central jet vetoes [87][88][89][90], wherein events that contain central (|η| < 2.5 − 3) jets with a transverse momentum above p Veto T are classified as background and rejected in searches and measurements.
While qualitatively sound, in reality, application of traditional jet vetoes with static values of p Veto T = 20 − 50 GeV do much to hurt signal selection efficiency in searches for heavy, new colorless particles [69,[98][99][100]. To quantify this statement, we plot in Fig. 14(a) the jet veto efficiency ε(p Veto T ) and uncertainty for the CC DY process, defined as the ratio:  over the mass range considered. This low acceptance can be attributed to the structure of initial-state radiation of gluons: after phase space integration, one finds that such contributions are of the form δσ/σ ∼ α s (µ r = p Veto T ) log 2 (Q/p Veto T ) 2 , for a hard scale Q ∼ m N , which indicates a propensity for heavier N to be accompanied by higher p T gluons. For example: the likelihood for a (relatively) soft gluon with momentump T to accompany a neutrino of massm N is essentially the same for a gluon of momentum 2p T to accompany a neutrino with mass 2m N . Hence, heavier N are accompanied by jets with p T that are more and more likely to exceed (and thus fail) the p Veto T threshold. (ii) In addition to this, one also sees that larger R jets uniformly fail the veto more readily than smaller R jets. This is directly due to larger R jets containing more constituents that are roughly traveling in parallel motion, and therefore contribute constructively to the jet's net p T . (iii) For this same reason, i.e., that larger R jets are more inclusive than smaller R jets, we observe a smaller scale dependence for large-R jets than small-R jets. However, this is also impacted by formally O(α 2 s ) terms in the NNLL resummation that scale as log(R), and therefore numerically vanish in the R → 1 limit, minimizing the jet radius dependence. Some unmatched scale dependence is also inserted by the matching of the NLO calculation to an NNLL resummation, and not a NLL resummation. For more technical details, see Ref. [96] and references therein. As noted in Ref. [100], this behavior is somewhat universal in that it is largely dependent on color structure and relevant mass scales, and hence can also be found when jet vetoes are applied to slepton/electroweakino pair production [98] as well as Sequential Standard Model W ′ /Z ′ boson production [100].
For higher collider energies, which is shown in Fig. 14 In the end, this shows that independent of collider, if one wishes to apply a jet veto with a fixed p Veto T threshold and retain sensible efficiencies, i.e., ε(p Veto T ) 80 − 95%, then one must restrict their searches to heavy neutrino mass scales of only a couple hundred GeV, or relax p Veto T to effectively act as if no veto is applied. While m N is an unknown quantity (assuming N exists), the premise of relaxing p Veto T with increasing m N to act as if no veto is applied is not exceptional. As discussed extensively in Secs. 4.3 and 4.4, the three leading charged leptons in the CC DY pp → N ℓ → 3ℓX process all possess momenta that scale with m N (see Eq. 4.37). At the reconstructed level, they do not appreciably change with variable collider energy and are distinguishable from continuum contributions (see Fig. 9). A similar argument for the leading charged leptons in the W γ → N ℓ fusion process also holds. Hence, in a genuine heavy neutrino event, the p T of the leading charged leptons scale proportionally with the neutrino mass and can act as a proxy for m N . Therefore, by setting p Veto T on an event-by-event basis, for example, to the leading charged lepton in an event, one would expect an increase in the jet veto efficiency for genuine heavy neutrino trilepton events.
Variable jet vetoes are not new, per se. In precision SM calculations [226][227][228], comparable choices for p Veto T have been made as a computational convenience to alleviate large veto logarithms of the form log 2 (Q/p Veto T ) 2 , rendering resummation unnecessary. Independently, the lepton-over-jet p T ratio (p ℓ T /p j T ) has also been used experimentally [229] at the LHC for particle identification to help distinguish leptons associated with hadronic / jet activity. However, as first reported in the companion letter of this report [69], we find that dynamic jet vetoes can be successfully used in a much broader class of experimental measurements and searches, notably searches for new, high-mass colorless particles.
We investigate the impact of a dynamic jet veto on the CC DY production of heavy neutrinos by first noting that the leading charged lepton in such events possesses a transverse momentum that scales as p ℓ 1 T ∼ m N /2. Hence, we set † the veto threshold to and plot in Fig. 14 Two additional features are observed: The first is a significant decrease in residual scale dependence, which spans from δε ∼ ±1% to ±5 for p Veto T = m N /2. Remarkably, the uncertainty exhibits little-to-no sensitivity to heavy neutrino mass scale, jet radius, and collider energy under this veto scheme choice. This is attributed simply to the fact that tying p Veto T with p ℓ T ∝ m N ∼ Q outright avoids the generation of large jet veto Sudakov logarithms; and for (m N /2) ≫ p Veto T = 30− 100 GeV, results in a much more inclusive cross section than in the static jet veto case. As a consequence of this reduction, we are able to observe a second feature, namely the appearance of a excess for the R = 0.1 curves, where one sees ε(p Veto T ) slightly surpassing unity for all √ s. Upon further investigation [69], we find that this originates from O(α 2 s ) terms in the NNLL(veto) resummation that scale as α 2 s log R. We confirm this by setting R = 0.01 and finding that ε(p Veto T ) ≈ 100 − 115% for m N 150 GeV, clearly illustrating a breakdown of the perturbative calculation and the need for a dedication resummation of log R terms as done in Refs. [92,201]. We stress that this does not indicate a breakdown of the dynamic jet veto scheme but only a breakdown of the NLO+NNLL(veto) computation in the R → 0 limit. The picture that emerges † It is not presently possible to set p Veto T on an event-by-event basis in the public NLO+NNLL(veto) code of Refs. [85,96], and hence we set the veto threshold here explicitly to p Veto T = mN /2. In all following parton shower-level discussions, p Veto T is set on an event-by-event basis to the appropriate quantity. from Fig. 14(b,d,f) may indeed hold for R ≪ 0.4, but further investigation is required and encouraged. For further details, see Ref. [69] and references therein.

VBF and jet vetoes
It is worth discussing how the W γ-fusion process, which possesses at least one energetic forward jet with p j VBF T M W /2 ∼ 40 GeV, is affected by a dynamic jet veto. An important conceptual point to stress is that jet vetoes do not require events to have zero hadronic activity. They select out only events with (central) jets possessing a transverse momentum above p Veto T but are inclusive with respect to jet activity outside this threshold. As a consequence, the VBF process will survive the veto as long as p Veto

τ h and jets vetoes
Another point worth examining is the treatment of hadronically decaying τ leptons (τ h ) in the presence of a jet veto. Experimentally, τ h are reconstructed first as jets before τtagging is applied [230,231]. Hence, one must justify whether a τ h can be excluded from the veto procedure, which requires demonstrating that τ h are objects totally independent QCD jets. As argued in Ref. [69], under the NWA, a final-state, on-shell τ is color-disconnected at a perturbative level from the rest of the hard pp collision, to all orders in α s . At a non-perturbative level, the τ → τ h ν decay occurs on a macroscopic distance of away from the pp interaction point. Here γ τ is the τ lepton's Lorentz factor, and Γ τ ∼ 2 × 10 −12 GeV is its total width. This is a much later transition than hadronization, which for a non-perturbative scale of Λ NP ∼ 1−2 GeV, occurs at a distance of d NP = 1−2 fm away from the pp collision. Consequently, τ leptons outlive primary hadronization of high-p T pp collisions, and the remaining coupling to the (pp)-system is through long-range, colorsinglet exchanges [169,[232][233][234], i.e., beyond twist-2 in the operator product expansion. Contributions of this kind are beyond the formal accuracy of the Collinear Factorization Theorem in Eq. 4.1, and hence can be consistently neglected, thereby demonstrating a decoupling of τ h from QCD jets in an event.

Dynamic Jet Vetoes and PDF Uncertainties
Owing to the complexity of jet vetoes, there are numerous sources of theoretical uncertainty. One such uncertainty is stemming from fitting PDF normalizations to data and its impact  on σ NLO+NNLL(veto) cross section predictions. For the CC DY pp → N ℓ + X signal process, we quantify this by considering the ratio of cross sections,  [235]; we also consider the NNLO NNPDF 3.0 set [236] as a representative "PDF4LHC Run II" recommendation [237]. Uncertainties are derived from PDF replicas in the appropriate statistical manner [176].
In Fig. 15 and as a function of heavy neutrino mass, we show the ratios (a) R NLO and (b) R NLO+NNLL(veto) (p Veto T = m N /2), with PDF uncertainties, at √ s = 14 TeV. In the upper (lower) plots, the MSTW 2008 reference curve is overlaid with the NNPDF 3.0 (NNPDF 3.1+LUXqed) curve. Notably, for all three PDF sets, we observe few differences, qualitatively and quantitatively, between the (a) no-veto and (b) jet-veto cases. For the mass range m N = 100 − 1000 GeV, we see that the uncertainty for all three PDFs span ±1% − ±3%. Qualitatively, we see that central value of the NNPDF 3.0 (NNPDF 3.1+LUXqed) PDF uniformly undershoots (overshoots) the reference PDF's central value but essentially always stays within its 1σ band. The slight exception is for the NNPDF 3.0 case at m N < 200 GeV, where the central value exceeds the 1σ lower band by ∼ 1%.
For the uncertainty at higher collider energies, we first note that the range of Bjorkenx considered here spans roughly x ∼ m N / √ s ≈ 0.01 − 0.07. This corresponds roughly to m N ≈ 300 − 1900 GeV at √ s = 27 TeV and m N = 1 − 7 TeV at √ s = 100 TeV, and hence covers the bulk of the investigated parameter space. Over these mass ranges, one expects a similar PDF uncertainties as those reported in Fig. 15 due to PDF scale invariance. For smaller x down to x ∼ (150 GeV/100 TeV) = 1.5 × 10 −3 , the NNPDF 3.1+LUXqed uncertainly is approximately unchanged [177]. For larger x up to x ∼ (20 TeV/100 TeV) = 0.2, the uncertainty grows larger but remains under ±5% [177].

Backgrounds and Jet Vetoes at NLO+PS
In light of the improved veto efficiency of the dynamic jet veto for the CC DY and VBF signal processes relative to the static veto, it is necessary to explore how background processes fare. For the inclusive pp → N ℓ + X → 3ℓ + X signal process, the leading SM processes fall into three categories: (i) top quarks, (ii) EW diboson and multi-boson production, and (iii) fake charged leptons, which we now discuss. Throughout Sec. 5.2, we note that to emulate minimal, analysis-level selection cuts all signal and background processes are evaluated at NLO+PS according to Sec. 3 and the following cuts are applied to jets and the three leading charged leptons after reconstruction:

Top Quark Background
Single and double top quark production processes at multi-TeV hadron colliders, e.g., (5.12) are major background for any multilepton measurement and search due to their large cross sections, intrinsic mass scales, and diversity of final states. Since BR(t → W b) ≈ 100%, multilepton search strategies for heavy neutrinos make use of b-tagging, and hence b-jet vetoing, to help suppress these backgrounds [34,74]. However, even for high-efficiency taggers, such as the CMS CSVv2 algorithm [238,239] used by Ref. [74], which possesses a tagging efficiency of ε b−tag ≈ 70 − 80%, at least 4 − 10% of top quark pairs survive the single-and double-tagging. Of the total number of top quarks, the fraction is actually larger when one takes into account that only jets within |η| 2.5, i.e., within tracker coverage, can be tagged. It has already been reported elsewhere [100] that relaxing the btagging requirement and simply vetoing central jets with p Veto T = 20 − 50 GeV, independent of flavor composition, can improve top quark rejection in searches for DY-like processes, even with an ideal efficiency of ε b−tag = 100%. Here, we improve upon this idea by allowing p Veto T to instead be set on an event-by-event basis to the p T of the leading charged lepton in our multilepton final state. On the other hand, the b-jet momenta, independent of being tagged, scale as Therefore, we see that the top quark backgrounds in multilepton searches for heavy N inherently fail the dynamic veto cut that requires p ℓ 1 T > p j 1 T . In fact, top quark processes fail regardless of choosing the leading charged lepton ℓ 1 , the subleading charged lepton ℓ 2 , or even the trailing charged lepton ℓ 3 .
To show this, we plot in Fig. 16 the (a,c,e) p T of the leading jet (p j 1 T ) and (b,d,f) leading lepton-to-leading jet p T ratio At √ s = 14 TeV, we see in that the CC DY signal and ttℓν background processes both exhibit their characteristic p j 1 T behavior at low-and high-p T respectively. Due to binning effects, the Sudakov shoulder in the signal process is unobservable. For ttℓν, the wide bump near p j 1 T ∼ 70−90 GeV corresponds to the anticipated p T for b-jets as given in Eq. 5.14. The higher value for the maximum is attributed to hard ISR and FSR that appears at O(α s ). Focusing now on the ratio plot, we see crucially a qualitative difference in the top quark and signal process: The signal process has an incredibly broad, continuum-like spectrum that appreciably starts at r ℓ 1 j 1 ∼ 0.25, possesses a very shallow peak at r ℓ 1 j 1 ∼ 1 (1.5) for m N = 150 (450) GeV, and readily spans rightwards by several units. The ttℓν background, on the other hand, peaks at a ratio of r ℓ 1 j 1 ∼ 0.5, with over 75% of events sitting below r ℓ 1 j 1 < 1. While not shown, we report that ratio for the subleading lepton r ℓ 2 j 1 = (p ℓ 2 T /p j 1 T ) features a slightly stronger separation between signal and background.
For higher collider energies, we observe encouraging behavior: While the p j 1 T distribution CC DY signal process broadens, with the lowest bin occupancy for m N = 150 (450 GeV) dropping from about 30% (22%) at √ s = 14 TeV, to 25% (17%) at 27 TeV, down to 20% (10%) at 100 TeV, the ttℓν rightward shift toward higher values of p j 1 T is more significant. This seen clearest in the r ℓ 1 j 1 ratios at 27 and 100 TeV, where about 35% and 40% of events, respectively, have a leading lepton-to-leading jet p T ratio less than r ℓ 1 j 1 < 0.5. This is in comparison to the roughly 25% of ttℓν events at 14 TeV. For the CC DY process, the migration of r ℓ 1 j 1 to values below unity is found to be only slight for m N = 150 GeV and mostly negligible for heavier neutrino masses.
At first, the dynamic jet veto choice of p Veto T = p ℓ 1 T does not appear to improve top quark discrimination over b-tagging central jets. We will see next, though, that in conjunction with additional information on charged lepton activity, it is sufficient.

EW Diboson Continuum and Resonant Multiboson Production
In multilepton searches for heavy neutrinos with masses at or above the EW scale, the production of two or more EW bosons, either non-resonantly or resonantly, represent the leading EW backgrounds at pp colliders [34,74]. However, despite their colorsinglet nature, they are significantly less immune to jet vetoes than the CC DY and W γ fusion signal processes. Due to the presence of radiation zeros in Born-level amplitudes, a large fraction of the inclusive diboson and triboson processes contain at least one high-p T jet [240][241][242][243][244][245]. For example: NLO studies of the inclusive pp → 3W + X process reveal that O(30%) of the cross section is comprised of the subprocess pp → 3W + 1j with p j T > 50 GeV [246,247]. As a result, NLO and NNLO corrections to EW diboson [248][249][250][251][252][253][254][255][256][257][258][259] and triboson [247,260,261] production uncover large (> +100%) increases in the total normalization of inclusive cross sections and a sizable increase in final-state jet multiplicity. This can be effectively observed in Fig. 17, where we plot the (normalized) number of jets with p T > 25 GeV at (a) √ s = 14, (c) 27, and (e) 100 TeV, for the CCDY signal process assuming m N = 150 (solid) and 450 GeV (dash) as well as the 3ℓν (dash-1 dot), W W W (dash-2 dots), and ttℓν (dash-3 dots) background processes. At √ s = 14 TeV, one sees that fewer than 40 − 45% of pp → 3ℓν and 3W events fall into the zero-jet bin, i.e., events with only jets possessing p j T < 25 GeV. This is comparable to the CC DY signal process at m N = 450 GeV but unlike the m N = 150 GeV process, which is closer to the 55-60% threshold. Conversely, the pp → 3ℓν and 3W processes more readily populate higher jet multiplicities, with 10% of events possessing at least three jets, whereas only 5% of CC DY events at m N = 150 GeV do. At higher √ s, one clearly sees that the lowest (highest) multiplicity bins deplete (grow) faster for the EW backgrounds than the CC DY samples. Concretely, for the signal process at m N = 150 (450) GeV, one observes that the fraction of events with zero jets above 25 GeV drop from about 60% (45%) at √ s = 14 TeV to roughly 50% (35%) at 27 TeV, to approximately 40% (25%) at 100 TeV. Similarly, for the ttℓν process, at least 10% of events have four or more jets at 14 TeV; at 27 TeV, at least 10% of events have five or more jets; and at 100 TeV, over 10% of events have at least six jets. Roughly speaking, at our level of modeling, about 5% of ttℓν events have eight or nine jets with p T > 25 GeV. One should caution in interpreting these results: For smaller jet radii, one expects the level of migration (and multiplicity) to be more dramatic. In addition, the accuracy here is only NLO+PS, indicating that the highest jet multiplicity bins are populated by the parton shower. To see how the EW background for the pp → N ℓ + X → 3ℓ + X process copes with a dynamic jet veto in light of this jet activity, we revisit Fig. 16, which additionally shows the representative pp → 3ℓν (dash-1-dot) and pp → 3W → 3ℓ3ν (dash-2-dots). As the EW processes above are driven by (qq)-annihilation, we see much the same qualitative behavior in the p j 1 T distribution as we do for the signal processes. Quantitatively, however, the larger jet activity of the multiboson channels causes both the 3ℓν and 3W processes, which have As a result, in the leading charged lepton-to-leading jet p T ratio r ℓ 1 j 1 = (p ℓ 1 T /p j 1 T ), we see that most of the EW events possess r ℓ 1 j 1 < 1. This is very unlike the CC DY process and in fact much more the like ttℓν process. As a function of collider energy, we observe a similar (though perhaps slightly milder) √ s dependence for r ℓ 1 j 1 as we do the ttℓν case, i.e., an increase separation of signal and background, and need not be discussed further.
An interesting consequence of the high jet activity in the EW pp → 3ℓν and 3W processes is the impact on the p T of the leading leptons. Unlike the CC DY and VBF signal processes, which feature low-p T ISR, the EW backgrounds feature high-p T ISR that recoil against the leading charged leptons. The impact of this can be seen in Fig. 17, where we show for (b) √ s = 14, (d) 27, and (f) 100 TeV, exclusive S T as defined in Eq. 4.35, for the several processes presently under discussion. For the EW processes, the p T of the charged leptons scale as p ℓ T ∼ M V /2, leading to a characteristic S T of This distribution is observed at √ s = 14 TeV, but with a particular broadness. We attribute the dispersion in S 2V T and S 3V T , which would otherwise be much narrower distributions since they result from from discrete transitions / decays, to the hard ISR. For reference, the broadness of the CC DY signal processes originate from the prompt charged lepton in the initial pp → N ℓ N scattering process, which leads to a continuum distribution for p ℓ N T (see Fig. 6). At higher collider energies, this broadening is worsened, which can be attributed to the associated increase in jet multiplicity. Now, with Fig. 17 in mind, imposing a dynamic jet veto of r ℓ 1 j 1 = (p ℓ 1 T /p j 1 T ) > 1 on the EW backgrounds would result in two outcomes: First is the removal of a sizable fraction of background events, as evident in Fig. 16. Second, is the reduction of hard ISR that presently kicks the multilepton systems and is responsible for broadening the S 3ℓν T and S 3W T distributions. This implies that the EW boson systems would largely be at rest and, crucially, that their decay products would possess Born-like kinematics more in line with Eq in conjunction with a dynamic jet veto would eradicate the EW (and top quark) backgrounds, while remaining resilient as a function of collider energy.

Fake Leptons
Non-prompt leptons from hadron decays and light-flavored QCD jets mis-identified as electrons or hadronically decaying τ leptons (τ h ), objects collectively known as "fake leptons," represent a non-negligible source of backgrounds in multilepton searches for heavy neutrinos at hadron colliders [47,59,69,[73][74][75][76]79]. Fake leptons satisfying selection criteria, however, must necessarily originate from high-p T hadronic activity. By color conservation, this implies [69] a high likelihood of additional hadronic activity, i.e., a second jet, with comparable p T elsewhere in the detector, and suggests that a jet veto can improve the rejection rate of fake lepton backgrounds. The reason for this is that high-p T jets and hadrons that give rise to fake leptons are seeded by high-p T partons that are charged under QCD. As the initial protons colliding form precisely a color-singlet state with zero, net transverse momentum, the color charge and transverse momentum of this parton must be balanced by one or more recoiling partons, which separately form hadrons / jets. For the specific case of hadron decays, e.g., B → Dℓν, this may simply be the spectator hadron, which is known to have similar momenta as the outgoing charged lepton [262,263]. Independent of this, rates for mis-identifying light jets as electrons and τ h are highest for the lowest p T jets in an event due to poorer resolution [173,174,264]. Hence, when taken together, the presence of a fake lepton implies not only the existence of one or more additional clusters of hadronic activity, but that these additional clusters have a slightly higher likelihood to possess a p T greater than the fake lepton itself. As a consequence, the use of a dynamic jet veto of r ℓ k j 1 = (p ℓ k T /p j 1 T ) > 1 for k = 1, . . . , 3, can improve the rejection rate of this class of background.

Dynamic Vetoes Beyond p T Ratios
At its heart, the dynamic jet veto we employ in Eq. 5.15, in conjunction with R = 1 jets, functions as a discriminant between leptonic and hadronic activities. The p T of leading leptons and jets, however, are not the sole measure of such activities. Other observables can also quantify this in complementary aspects due to their variable level of inclusiveness / exclusivity. For example: For leptons, there is exclusive S T , as defined in Eq. 4.35. For hadronic activity, there is inclusive H T , which is defined as the scalar sum of hadronic p T , Therefore, in this section, we briefly explore if other dynamic selection criteria can also fulfill our main intent of improving sensitivity to multi-lepton searches for heavy N . As in Sec. 5.2, signal and background processes are evaluated at NLO+PS according to Sec. 3 and the nominal fiducial and kinematic cuts of Eq. 5.9 are applied to jets and the three leading charged leptons after reconstruction. The premise for considering alternative measures follows from the observation in Sec. 4.4 that there exists a class of leptonic observables, not just p ℓ k T for k = 1, . . . , 3, whose distributions are largely insensitive to varying collider energies. Such a property, though, does not hold in general, as seen throughout Secs. 4.4 and 5.2. For example, in Fig. 17, for all processes there is an uptick in the jet multiplicity with increasing √ s due, in part, to the opening of phase space. A consequence of increasing jet multiplicity is the increase in hadronic energy carried away from the hard scattering process to the detector experiments.
To quantify this, we show in Fig. 18 the normalized inclusive H T distribution, as defined in Eq. 5.20, at (a) √ s = 14, (c) 27, and (e) 100 TeV, for the CCDY signal process assuming m N = 150 (solid) and 450 GeV (dash) as well as the 3ℓν (dash-1 dot), W W W (dash-2 dots), and ttℓν (dash-3 dots) background processes. Interestingly, while one sees a considerable broadening of the H T for both signal and background processes, the background processes broaden at a slightly faster rate and the broadening is most pronounced for the 3ℓν and ttℓν. Moreover, for increasing √ s, a rightward lurch to larger H T can be observed in the ttℓν. The acute sensitivity of H T to collider energy suggests that, like the ratio p ℓ 1 T /p j 1 T , the ratio p ℓ 1 T /H T may serve as criterion to base a veto on hadronic activity. On the other hand, the robustness of S T over varying collider energy for the signal process suggests that the ratio S T /H T may too serve as a comparable, if not better, discriminant of leptonic and hadronic activity. To investigate this, in Fig. 18(b,d,f), we at last show the ratio S T /H T for √ s = 14, 27, and 100 TeV, respectively. Remarkably, as √ s increases, the background processes grow significantly more narrow than the CC DY signal processes and shifts leftward to smaller S T /H T , suggesting a potentially powerful means to reject backgrounds at future colliders that inherits the single-scale properties of the p ℓ 1 T /p j 1 T discriminant. Further optimization of jet veto selection criteria, such the relative gain (or lack thereof) of choosing S T /H T > 1 as a veto criterion rather than p ℓ 1 T /p j 1 T > 1, are beyond the scope of this work and are deferred to future studies. Similarly, investigations into whether the precise requirement that r ℓ 1 j 1 = p ℓ 1 T /p j 1 T > 1 is better or worse than choosing, for example, r ℓj > 0.5 or r ℓj > 2, are strongly encouraged, particularly in the context of a multivariate analysis. Nevertheless, for the remainder of this study, we choose as our dynamic jet veto threshold r ℓ 1 j 1 > 1.   Figure 19.

Dynamical Jet Vetoes at Leading Logarithmic Accuracy
As discussed in Sec. 5.1, a principle consequence of choosing a dynamical p Veto T in lieu of a static veto is a reduced dependence of jet veto cross sections on IR and UV cutoff scales, e.g., µ f , µ r , µ s . This is clear at NLO+NNLL(veto) in Fig. 14, where scale dependencies reduce from the 1-to-15% level to the sub-percent level for DY production of N . This subsequently raises two questions: (i) Does such a reduction also holds for other color structures and scattering topologies (and hence different radiation patterns)? (ii) How large is the scale dependence at lower logarithmic accuracies, particularly at LL, which is universally calculable with parton showers? Finding a universal reduction in scale uncertainties for jet veto cross sections of processes with alternative structures and at formally lower precision would facilitate a broad application of dynamic jet vetoes in experimental searches. Such a broad, systematic investigation, however, is beyond the present study but encouraged. Instead, we take a more limited but highly illustrative first step.
For the DY and VBF signal processes at representative heavy N masses and for representative background processes (pp → 4ℓ, 3ℓν, ttℓℓ, and ttℓν), we report the scale dependence of jet veto cross sections and efficiencies at √ s = 14 TeV, up to NLO+PS(LL), respectively, in Tables 2, 3, and 4, and graphically in Fig. 19. Specifically, we determine the factorization (µ f ) and renormalization (µ r ) scale dependence in cross sections when jointly scanned over a three-point scale variation † , at LO+PS and NLO+PS. As a normalization, we define the a generic particle-level, fiducial cross section (4th column) obtained from applying the same requirements as given in Eq. 5.9 on jets and the three leading charged leptons in an event. The total inclusive rate, i.e., without kinematic selection criteria beyond necessary generator-level cuts, is also reported (third column) for completeness. In addition to Eq. 5.9, we consider the case of a dynamical veto with p Veto   Table 2. For the CC DY pp → N ℓ + X → 3ℓν + X signal process, with e − µ/no-τ mixing as given in Eq. column) as well as the case of a static veto with p Veto T = 30 GeV (sixth column). For alternative √ s, we conjecture that one should observe comparable changes in scale dependence as observed in Fig. 14 and Ref. [100]. To quantify the impact of QCD corrections, we report the NLO+PS K-factor at each cut iteration. The statistical confidence correspond to 500k generated events, except for the ttℓℓ samples, which correspond to 1M events.
For the DY signal process (Table 2), a number of qualitative features can be discerned: Foremost is that the total inclusive, fiducial, and both jet veto scale uncertainties at NLO+PS are all comparable in size, and span at most ±2%; at LO+PS, the uncertain-  Table 3. Same as Table 2 but for the VBF pp → N ℓj + X → 3ℓνj + X signal process.
ties are 2-to-10× larger, but do not exceed ±7%. For the dynamical veto, efficiencies at NLO+PS span approximately 90-98% for the dynamical veto and 45-65% for the static veto. For increasing m N , veto efficiencies increase (decrease) under the dynamical (static) veto. Importantly, this is qualitatively and quantitatively comparable to both the central values and uncertainties observed at NLO+NNLL(veto), in Fig. 14. This suggests that NLO+PS provides a good description of dynamical and static jet vetoes, and sufficient for discovery purposes. For static vetoes, this is consistent with the findings of Refs. [97,100]. For dynamical vetoes, this is the first such quantitative comparison between parton showers and higher logarithmic jet veto resummations. In addition, the spectrum of K-factors indicate the relative contribution of real and virtual radiation in heavy lepton production via the DY mechanism: At larger (m N / √ s), the inclusive, fiducial, and both veto K-factors are comparable, indicating that QCD corrections are dominantly virtual. At smaller (m N / √ s), more exclusive / less inclusive final states possess smaller K-factors, indicating the increasing presence of real radiation. Owing to the uniformity of the NLO+PS K-factors for the dynamical veto, it appears that modeling the DY signal processes at LO+PS using the scale scheme of Eq. 3.16 and normalizing by a multiplicative K-factor of K ≈ 1.1 is a reasonable procedure. Such a prescription does not hold for the static veto. For the W γ signal process (Table 3), the overall behavior is comparable to the DY channel. First we observe that the inclusive, fiducial, and dynamical veto rates exhibit a common O (5) For the static veto, the situation is more dismal, with efficiencies of O(10)% due to the p T of the leading jet scaling as p j 1 T ∼ M W /2, which is above the threshold. While the precise choice of p Veto T = 30 GeV and the |η|-window are somewhat arbitrary, the principle at hand remains: for the mass scales under consideration, VBF topologies are not robust against arbitrary static jet vetoes due to the presence of forward, high-p T jets. We conclude that modeling the signal process at LO+PS with a QCD K-factor of 1 − 1.05 for a static veto and 1.05−1.1 for dynamic veto is a sufficient description of the VBF processes at NLO+PS.
Turning to our representative backgrounds (Table 4), the situation is more complex. Beginning with the pp → 3ℓν process (first row), one notices first the large difference between the NLO+PS and LO+PS cross sections for the inclusive, fiducial, and dynamical veto rates. The K-factors span K ≈ 1.  [248-251, 258, 259]. The huge increase is due to a radiation zero [240][241][242][243][244][245] in the qq ′ → W γ ( * ) process, which suppresses the Born rate, and is broken in the qq ′ → W γ ( * ) g and qq ′ → W γ ( * ) gg sub-processes. For the static veto, the K-factor is K ≈ 1.2, indicating the modest size of virtual and unresolved real radiation. For the dynamical and static vetoes, the efficiencies are about 75% and 50%, respectively, showing the utility of jet vetoes even on EW background processes. Due to a much weaker radiation zero, the pp → 4ℓ process (second row) exhibits a smaller K-factor but much the same scale uncertainties at NLO+PS and LO+PS. Therefore, aside from a slightly larger veto efficiency of ε ≈ 90% and 65% for the dynamical and static veto, respectively, the channel needs not to be discussed further.
Qualitatively, the top quark channels pp → ttℓν (third row) and ttℓℓ (fourth row) exhibit a very different scale dependencies and veto efficiencies than the 3ℓν and 4ℓ processes. The pp → ttℓν (ttℓℓ) process reveals K-factors of K ≈ 1.45 − 1.6 (1.45 − 1.75) across all rates, with O(20 − 30)% uncertainties at LO+PS that reduce to O(10)% at NLO+PS. At NLO+PS, the dynamic (static) veto cross section possesses a scale dependence a bit lower (larger) than this average; the difference between the two is about 1.7 − 2×. Notably, in the absence of other selection criteria, the dynamic and static veto efficiencies differ radically, with about 30% and < 1%, respectively. The very low static veto efficiency is due almost entirely to the existence of two high-p T b quarks in the decays of the tt pair. The much higher efficiency (but still low in absolute terms) observed for the dynamic veto is due to the prompt (ℓν) system is recoiling against a larger multi-body system. Roughly, boosting this prompt charged lepton by an additional 10 GeV translates to only a 5 GeV recoil to each b quark, and therefore permits the low-p b T tail of events to pass the veto. This suggests, however, that alternative dynamical choices for p Veto T , such at p Veto T , could further lower the dynamical veto efficiency. (We have verified this but find the improvement unnecessary due to other cuts applied.) Numerically, the ttℓℓ process contains larger K-factors after cuts than the ttℓν. This is attributed to the presences of contributions like gg → ttℓℓ, which are present at the Born-level, and receive larger virtual corrections than the ttℓν process, which is strictly initiated by (qq ′ )-annihilation at LO.
Collectively, the results of this discussion are summarized and illustrated in Fig. 19 for both the (a) static and (b) dynamic jet vetoes.

Observability of Heavy Neutrinos at Hadron Colliders
In Secs. 4.4 and 5, we established the existence of measures for lepton and hadronic activities that remained robust to changes in collider energies. We also showed that, when used together, namely a dynamic jet veto in conjunction with exclusive S T , the two could greatly improve signal acceptance and background rejection. We now turn to quantifying the impact of such a dynamic jet veto on the discovery potential of trilepton searches for heavy neutrino. The remainder of this section continues as follows: In Sec. 6.1, we describe our collider detector modeling, including object definition requirements and tagging/misidentification rates. We continue in Sec. 6.2 with defining our proposed dynamic jet veto analysis and benchmark trilepton analysis, which is based closely on the √ s = 13 TeV CMS analysis reported in Ref. [74]. Finally, in Sec. 6.3, we report our results for search prospects for Dirac and Majorana neutrinos at the √ s = 14 TeV, as well as projections for a 27 TeV HE-LHC and a hypothetical 100 TeV VLHC.

Hadron Collider Detector Modeling
In this section, we describe our modeling of a generic LHC detector experiment. The fiducial volume and segmentation are based on the ATLAS [265][266][267] and CMS [268][269][270] detector experiments. Detector response modeling, particle identification (PID) tagging and mistagging rates, as well as kinematic and isolation acceptances/thresholds, are based on published physics analyses, dedicated calibration and detector performance studies, where publicly available, and published trigger menus. All of which we now summarize.

Global Fiducial Volume
Our analysis starts from the particle-level event samples described in Sec. 3. Particle-level objects are stable on a detector's length scale and assumed, momentarily, to be reconstructed with 100% efficiency if they fall within a subdetector's fiducial volume. All objects with ultra forward rapidities, i.e., |y| > 5, are ignored outright; electrons, muons, hadronically decayed tau leptons (τ h ), and photons that fall outside the electromagnetic calorimeter (ECAL) coverage, i.e., |y| > 3, are relabeled categorically as light jets.

Detector Response
We next smear the momenta of these objects according to their (potentially new) PID. For all objects we employ a Gaussian smearing profile. Explicitly, this means that for some kinematic observableÔ, e.g.,Ô = p T or E, its value is perturbed randomly according to a normal distribution centered on the original value ofÔ with a spread of σÔ, i.e., In this language, the resolution (δ) ofÔ is the (dimensionless) quantity, δ = σÔ/Ô. The 4momentum is then recalculated assuming that the direction of a relativistic particle always remains unchanged. This is because the direction of an infinitely energetic stable object can be measured in hermetic detectors with high certainty, unlike its energy. Momentum reconstruction for electrons, photons, jets, and hadronic taus is determined largely through calorimetry. Hence, we smear their momenta via shifts in their energies. For electrons and photons, the energy smearing is parameterized by [138,271,272] σ Here b e E = 2% for E T < 500 GeV at all η. For larger E T , b e E = 4 (6)% for objects inside (outside) the central barrel, which extends to |η| < 1.444. For simplicity, electrons and photons are treated identically and so we set b γ E = b e E . For jets, we combine the energy smearing adopted in the 13 TeV pp → tt+nj analysis of Ref. [273], which exploits energy resolution measurements for R = 0.5 − 0.7 jets [273,274], with p T smearing based on a dedicated calibration of R = 1.0 jets [275]. Jet energies and p T are varied independently so that changes in a jet's momentum are translated into a shift in the jet's mass, again leaving the 3-momentum direction unmodified. Like electrons and photons, the jet energy smearing is given by where the coefficient is b j E = 3% (5%) for central (forward) jets with rapidities satisfying |y| < 3 (|y| > 3) [274]. From Ref. [275], one can extract the jet p T smearing function For jets with p j T ∈ [250 GeV, 1.5 TeV], Ref. [275] reports that the power-law expression and the values of the coefficients are nearly uniform over the most central region of the ATLAS detector's barrel. We therefore extrapolate this to all jets with |y| < 3. Outside this p T window, we assume a linear function with coefficients set by the boundaries of Eq. (6.4), where b j p T ≈ 7.6% (3.2%) for p j T < 250 GeV (p j T > 1.5 TeV). For forward jets with |y| > 3, we use a blanket, y-independent coefficient of b j p T ≈ 7.6%. As τ h 's are experimentally a special class of jets, we apply the same smearing protocol to them as we do to QCD jets.
The momenta of muons are determined via track curvature. Hence, their smearing is modeled correspondingly through shifts in p T , with deviations (σ p µ T ) parameterized by.

PID Tagging and Mistagging
Summarily, PID tagging (ǫ Tag. ) and mistagging (ǫ Mis−Tag. ) probabilities are estimated from published ATLAS and CMS physics analyses and detector performance studies. For a given object, PID tagging/mistagging is implemented simply by comparing a randomly generated number (with TRandom3::Uniform()) to the probability, which may depend on p T . We first check if b-jets survive the tagging cull using the DeepCSV (loose) efficiencies as a function of jet p T , as reported in Ref. [239]; those that do not survive are henceforth identified as light jets. Next, τ h are checked, using the p T -dependent eτ h -Era-2017-F tagging efficiencies of Ref. [264]; those τ h that are not tagged are identified as light jets. τ htagging rates span (roughly) ǫ τ h →τ h ∼ 35−95% for p T 30 GeV [264]. After, genuine lightflavor QCD jets (and charged leptons and photons that are outside the ECAL coverage) are checked if they are misidentified first as (a) b-jets according to Ref. [239], the rates of which span ǫ j→b ∼ 1 − 3% for p j T 20 GeV; then as (b) τ h assuming the p T -dependent Tight MVA Discriminant rates of Ref. [174], which span ǫ j→τ h ∼ 0.01 − 1% for p j T 20 GeV. The electric charge of jets misidentified as τ h is assigned with equal probability. Those light-flavored jets that escape being mistagged are grouped with genuine b-jets and τ h that were misidentified as light jets.
We take into account the possibility of a light QCD jet being misidentified as an electron/positron, which is a main component of the "fake lepton" background in heavy neutrino trilepton searches [34,74]. For those background processes where such an occurrence is potentially important (see Sec. 3.3), on an event-by-event basis, we randomly choose a jet from the event (using TRandom->Integer()) and relabel it either an electron or positron (with equal probability). The likelihood for the event itself is then re-weighted uniformly by a factor of ǫ j→e = 7.2 × 10 −5 [173], independent of jet kinematics. Throughout the analysis, electron/positron charge mis-measurement is neglected.

Isolation and Analysis-Object Definitions
Stable electrons and muons (ℓ ± ) are considered hadronically isolated when the sum of the total hadronic E T within a distance of R max centered on the the lepton candidate is less than a fraction ε ℓ Had. Iso. of its E T , i.e., Had. Iso.
for ∆R ℓX < R max . (6.7) Following Ref. [74], we use a loose isolation criterion and set ε ℓ Had. Iso. = 30%, R max = 0.3. Reconstructed jets, and hence τ h , are inherently isolated from other pockets of hadronic activity since they are built up from sequential jet clustering algorithms [276]. Furthermore, two leptons (ℓ 1 , ℓ 2 ) are considered leptonically isolated if they satisfy the separation, Analysis-quality charged leptons and jets are subsequently defined as isolated objects satisfying the following fiducial and kinematic criteria: We stress that in hadron collisions, not all isolated objects satisfy the above requirements.
In particular, stray QED emission off electrically charged leptons and partons can give rise Table 5. Signal categories for trilepton signal processes mediated by a heavy Dirac neutrino (N ), underlying mixing hypotheses, whether the signal process is charged lepton flavor-conserving (cLFC) or -violating (cLFV). Here ℓ X ∈ {e, µ, τ h }, l i , l j ∈ {e, µ}, and no ± indicates that both lepton charges are permitted.
to central, soft electron and muon pairs; relatively soft muon pairs and τ h can also originate from decays of hadronic resonances; and hard, non-diffractive hadron collisions inherently give rise to forward, low-p T jets due to color conservation and confinement. Independent of the multiplicity of analysis-quality objects, and as a proxy for the net impact of energetic, final-state, light neutrinos, we define the transverse momentumimbalance vector ( p T ) and its magnitude (MET) by the following: The summation here is over all visible objects within the fiducial volume of the detector. The exception to this are reconstructed jets with p T < 1 GeV.

Signal Process and Collider Signature Definitions
We now describe our proposed dynamic jet veto and benchmark analyses at √ s = 14 TeV.
Dynamic Jet Veto Trilepton Analysis at the √ s = 14 TeV LHC As our underlying signal process, we consider the inclusive production of a single heavy neutrino N and a charged lepton ℓ N through the CCDY and W γ VBF processes, with N decaying through a SM charged current to a fully leptonic final state, given by In accordance to the benchmark active-sterile flavor mixing scenarios, given in Eqs. 3.17 and 3.18, we define several flavor permutations for our signal process, which we categorize by flavor hypothesis and summarize in Tab. 5. For a particular flavor signal category, we define our collider signature as precisely three analysis-quality charged leptons and MET, pp → ℓ 1 ℓ 2 ℓ 3 + MET + X, for ℓ k ∈ {e, µ, τ h }. (6.13) Here the charged leptons ℓ k (as well as all other objects) are ordered according to their p T , with p k T > p k+1 T , and, as above, X denotes an arbitrary number (including zero) of central, high-p T jets. To minimize contamination from multi-EW boson processes, we reject events with four or more analysis-quality charged leptons, but remain inclusive with respect to charged leptons that do not meet the analysis object definitions in Sec. 6.1.
We remove leptons from the Z pole and the low-mass SM Drell-Yan spectrum by requiring the following invariant mass cuts: (6.14) Dilepton invariant masses are built from all possible (ℓ i ℓ j ) permutations, independent of flavor or electric charge, to suppress electric charge mis-measurement and fake leptons. The trilepton invariant mass suppresses contributions from rare, but nonzero Z → 4ℓ decays.
To curb EW, top quark and fake backgrounds, we impose a dynamical central jet veto, namely the so-called "safe jet veto" [69], and set on an event-by-event basis More precisely, we require that all analysis-qualitatively jets possess p T less than the p T of the event's leading charged lepton. Events with any number of analysis-quality jets possessing a p T greater than the event's leading charged lepton are cut. We reiterate that the jet veto does not eliminate all jet activity in the event. We remain inclusive with respect to soft and forward hadronic activity: events may contain an arbitrary number of jets with p T < max(25 GeV, p ℓ 1 T ) and/or |y| > 4.5. (We briefly report that a slight improvement in the S/B ratio was observed when setting p Veto T = p ℓ 2 T ; the effect was less pronounced for the trailing charged lepton. We encourage further investigation into the matter.) To further repress EW and top quark production, we demand for ℓ 1 , . . . , ℓ 3 that S T > 125 GeV. (6.16) This requirement severely impacts the survival of continuum 3ℓν and 4ℓ processes since such cross section scale inversely with S T , with σ(pp → nℓ + X) ∼ 1/M 2 nℓ ∼ 1/S 2 T . As a proxy to the mass of N , we build the multi-body transverse mass variable (M M T ),  The dependence of the mass-window on m hypothesis N reflects two realities: (i) The total width of TeV-scale N can be numerically large (though still perturbative), since Γ N ∼ G F m 3 N |V | 2 , which intrinsically broadens the true value ofM T . As we consider m N ≥ 150 GeV, the mass window is never smaller than +15 GeV −20 GeV . (ii) Experimental resolution and the presence of MET also feed into the broadening of reconstructedM T . The asymmetrical requirements reflects the asymmetric nature of multi-body cluster mass variables; see Fig. 13. The precisely boundaries of −15% and +10%, however, are somewhat arbitrary and can be tuned. In practice, application of the above selection cut is no different from mass hypothesis-dependent diphoton or 4-charged lepton invariant mass requirements used in searches for the SM-like Higgs boson [277,278].
We summarize the above selection cuts in the top two rows of Table 6.
"Standard" Trilepton Analysis at the LHC Searches for heavy neutrinos in multi-charged lepton final states have long been wellmotivated since they compliment dilepton, one-lepton, and zero-lepton search channels. Subsequently, strategies premised and centered on the existence of high-p T charged leptons, such as those as proposed in Refs. [34,35,46,77], are now standardized; see, for example, the LHC searches of Ref. [74]. The analysis we propose qualitatively differs from these studies in several respects: (i) The primary difference is our use of a flavor-independent, dynamical jet veto; aforementioned studies consider at most a veto on b-tagged jets. (ii) Our requirements on individual charged leptons are relatively lax and uniform; previous studies tend to employ different, but nonetheless stringent, kinematic criteria for leading and non-leading charged leptons. (iii) We set stringent kinematic requirements on global leptonic activity. In a sense, our analysis takes advantage of and discriminates against the differences in global hadronic and leptonic activities of signal and background processes in hard pp collisions. This analogy is made firmer by our use of R = 1 jets.
It is therefore useful to quantify how our proposed analysis can improve (if at all) the anticipated LHC sensitivity. To this extent, we also consider an alternative "standard analysis" based closely on the CMS collaboration's trilepton heavy neutrino search at Run II of the LHC [74]. We define this benchmark analysis assuming the same lepton flavor collider signature of Eq. (6.13). After imposing the same leptonic invariant mass cuts of Eq. (6.14), we impose the additional kinematic criteria on an event's three charged leptons: We refer readers to Ref. [74] for the justification of these cuts. Events with at least one analysis-quality, b-tagged jet are vetoed. The "standard analysis" selection cuts are summarized in bottom row Table 6. Remarkably, under the same active-sterile mixing hypothesis as Ref. [74] and the same integrated luminosity at √ s = 14 TeV as reportedly used for 13 TeV (L ≈ 36 fb −1 ), we find good agreement with their expected sensitivity. Our implementation is slightly less sensitive than Ref. [74]. This serves as a highly nontrivial check of our signal and background modeling, including "fake" leptons, and our detector modeling.

Results: Sensitivity at the LHC and Beyond
In this section, we report the main results of our study, namely, the anticipated sensitivity to heavy neutrinos via the trilepton signature at current and future hadron colliders. Our findings are organized according to the following: In Sec. 6.3.1, we present the sensitivity to heavy Dirac neutrinos at the 14 TeV LHC, under various active-sterile mixing hypotheses, using our proposed analysis based on the dynamical jet veto as well as the standard, benchmark analysis. In Sec. 6.3.2, we repeat the exercise but for heavy Majorana neutrinos. Finally, in Sec. 6.3.3, we show the sensitivity to heavy Dirac neutrinos at the LHC's proposed 27 TeV upgrade and a hypothetical 100 TeV successor collider. In all cases, we report the sensitivity assuming Gaussian statistics. That is, for N s (b) signal (background) events expected with an integrated luminosity of L and a cross section of σ All Cuts after all selection cuts are applied, the signal significance (S) is quantified by . (6.20) A background systematic factor of δ b = 0.1 is applied to account for mismodeling of SM background processes, e.g., missing higher order QCD corrections and subleading processes, and detector effects, e.g., non-perfect electron and muon identification. For given L and σ All Cuts b , a particular mass and mixing hypothesis (m N , |V ℓN | 2 ) can be falsified at the (approximately) 95% confidence level (CL) if the corresponding signal rate σ All Cuts s results in a significance of S > 2. Fixing instead S 95 = 2, we can invert this inequality into an 95% CL upper bound on σ All Cuts s at a given L, which we label σ 95 , Using the relation between the bare cross section and mixing S ℓℓ ′ given in Eq. 4.6, the upper bound σ 95 can then be translated into an upper bound on S ℓℓ ′ at the 95% CL: where σ 0 is the bare cross section as derived from σ All Cuts s and any of the mixing hypotheses S Hypo. Table 5. We report our results in terms of the 95% sensitivity to S ℓℓ ′ , or equivalently |V ℓN | 2 given our mixing hypotheses, as a function of heavy neutrino mass m N . Throughout this section, we report the 95% CL sensitivity to active-sterile mixing |V ℓ4 | as a function of heavy neutrino mass N using the proposed dynamic jet veto trilepton analysis (solid) and the benchmark trilepton analysis (dash) for the signal categories as defined in Table 5. We use the mixing hypotheses as given in Eq. 3.18, which assumes N couples to only two charged leptons with equal strength, and Eq. 3.17, which assumes N couples to only one charged lepton. The two mixing scenarios are referred to as the charged lepton flavor violation (cLFV) scenario and the charged lepton flavor conservation (cLFC) scenario, respectively. We assume the nominal LHC and HL-LHC luminosity benchmarks of L = 300 fb −1 (darker, upper curves) and 3 ab −1 (lighter, lower curves) at √ s = 14 TeV.
Also shown for references are limits on |V e4 | and |V µ4 | from direct searches for the trilepton process by the CMS experiment at √ s = 13 TeV using L ≈ 36 fb −1 of data [74].
In Fig. 20, we plot the anticipated sensitivity to the trilepton process for the cLFV scenarios: (a) EMU-I, (b) EMU-II, (c) ETAU-I, (d) ETAU-II, (e) MUTAU-I, (f) MUTAU-II. We observe several common features: (i) For the lowest masses considered (m N = 150−300 GeV), we find that the dynamic jet veto analysis improves the sensitivity only marginally compared to present strategies. This is attributed, in part, to the stringent S T cut in Eq. 6.16, which greatly overlaps with the signal region for the lightest N , and to our preferred choice of transverse mass variable mass (see discussion near Eq. 4.40). (ii) For the highest masses (m N = 1 − 2 TeV), we find that that the new analysis improves sensitivity by roughly 7 − 15× at L = 300 fb −1 and well over 10× at L = 3 ab −1 . (iii) This increase in improvement is so large that for much of the mass range considered, the dynamic veto analysis with L = 300 fb −1 performs better than the standard analysis at 3 ab −1 . (iv) For most cases, we see that mixing as low as |V ℓ4 | 2 ∼ 3 − 6 × 10 −4 can be probed for m N 300 GeV and as low as ∼ 5 × 10 −3 for m N ∼ 1 TeV. The exceptions are (d) ETAU-II and (f) MUTAU-II, which include fewer sub-channels than other signatures and also suffers from double and triple τ h -tagging. For these channels, the sensitivity is uniformly reduced by about one order of magnitude. The similarity in reach across the various flavor scenarios follows from the near identical nature of the collider signatures themselves (see Table 5), meaning that differences are due to particle identification. Furthermore, for m N 300 GeV, the improvement of the trilepton analysis now makes it competitive in sensitivity to hadron collider searches for LNV mediated by heavy Majorana neutrinos [34][35][36][37][38], which possess considerably fewer backgrounds than the trilepton signature.  Table 5. Also shown are limits from direct LHC searches [74].  Table 5.
In Fig. 21, we plot the anticipated sensitivity to the trilepton process for the cLFC scenarios: (a) EE, (b) MUMU, (c) TAUTAU-I, and (d) TAUTAU-II. Qualitatively, the relative improvement of the dynamic jet veto analysis over the standard analysis is largely the same. An interesting exception to this is (d) the TAUTAU-II signature, where the performance of the standard analysis is noticeably better for m N 300 GeV. Here, we believe that a dynamic veto of p Veto T = p ℓ 1 T actually does more harm than good since all three charged leptons carry significant less momentum than the scaling behavior summarized in Eq. 4.37. For τ → τ h + ν decays, one expects τ h to carry only about half the momentum of the parent τ lepton; likewise, for τ → e/µ + 2ν, the e/µ carries only about a third of the original momentum. This implies that the p Veto T threshold is much lower in reality, and reduces signal acceptance, as see in Fig. 14. Quantitatively, we see a uniform reduction in sensitivity by nearly a factor of 2× with respect to the cLFV cases. This is directly tied to a smaller pp → N ℓX production cross section for N coupling to only one charged lepton, the difference being σ Tot. = σ(N ℓ 1 X) as oppose to σ Tot. = σ(N ℓ 1 X) + σ(N ℓ 2 X). In Fig. 22, we plot the anticipated sensitivity to the trilepton process mediated by a single Majorana neutrino for the cLFV scenarios: (a) EMU-I, (b) EMU-II, (c) ETAU-I, (d) ETAU-II, (e) MUTAU-I, (f) MUTAU-II. In all cases, we observe very comparable sensitivity for the Majorana neutrino scenario as we do for the Dirac neutrino scenario. For the benchmark trilepton analysis, the two sets of results are nearly indistinguishable, which follows from the analysis being largely independent of the heavy neutrino's Majorana character. For the dynamic veto analysis, the Majorana case features a slightly worse sensitivity, which we attribute to the increased likelihood of building the incorrect M M T .
In Fig. 23 we plot the anticipated sensitivity to the trilepton process for the cLFC scenarios: (a) EE, (b) MUMU, (c) TAUTAU-I, and (d) TAUTAU-II. Again, little difference between the Dirac and Majorana cases are observed and need not be discussed further.  Table 5. Also shown are limits from direct LHC searches [74], indirect global constraints [124], and perturbativity.  Table 5.
weaker in relative reach due to stronger limits on µ − τ mixing.
In Fig. 25, we show the results for the charged lepton flavor conserving signal categories (a) EE, (b) MUMU, (c) TAUTAU-I, and (d) TAUTAU-II. As seen for the √ s = 14 TeV in Sec. 6.3.1, the reach is weaker for the cLFC channels due to the smaller production cross sections but still of the same order of magnitude.

Summary and Outlook
Heavy neutrinos (N ) represents one of the best-motivated (though far from only) solutions to how active neutrinos are so much lighter than those elementary fermions whose masses have been confirmed to originate via the SM Higgs mechanism [279][280][281]. If they exist, however, their masses and the degree to which they couple to SM particles are far from clear. Hence, one must take a broad approach in searching for heavy neutrinos, particularly at collider experiments, which can directly probe N with EW-and TeV-scale masses.
To this extent, we have systematically reassessed the discovery potential of heavy neutrinos with masses m N ≥ 150 GeV, decaying to purely leptonic final states at the √ s = 14 TeV LHC, as given, for example, in Fig. 2. The impetus for this examination is two-fold: (i) Present LHC searches for heavy neutrinos [73][74][75][76] follow the seminal hadron collider analyses of Refs. [34,35,46,77], which justifiably argue that search strategies should rely on the existence of high-p T charged leptons when N are produced solely through the charged current Drell-Yan (CC DY) mechanism. It is now known, though, that new production mechanisms can compete or outright dominate over the CC DY mechanism; new tools have been created that can more reliably model heavy neutrinos in hadron collisions; and the understanding of jets has so significantly improved that jet observables can be used reliably as discriminates in measurements and searches. (ii) In addition to this, the community is presently assessing possible successors to the LHC program [27,[29][30][31]33]. Understandably, the analyses prescribed in Refs. [34,35,77] were not designed to be robust against increasing collider energy. It is unclear, therefore, if future projections based on current search strategies can reliably estimate the actual discovery potential of a √ s = 27 HE-LHC or a √ s = 100 TeV VLHC. We report investigating the matter and propose a new search strategy that both improves current sensitivity and remains robust at future colliders. To reach this result, we have considered the following: In Sec. 4 we have explored and compared the resonant production of heavy neutrinos N through a variety of EW mechanisms for √ s = 14 − 100 TeV and m N ≥ 150 GeV. At 14 TeV, the CC DY and W γ fusion (VBF) channels are dominant; at 27-50 TeV, gluon fusion (GF) emerges as an important and complementary channel for m N ≈ 300 − 1000 GeV; and outright dominates for m N 2 TeV at 100 TeV, beyond which VBF is the leading production vehicle. We stress that few, if any, studies exist demonstrating the discovery potential of such heavy N through the GF channel at hadron colliders. Normalized distributions were then presented for a number of leptonic observables across √ s = 14 − 100 TeV at LO and NLO+PS precision. This revealed the existence of a class of observables, namely exclusive quantities built from the transverse momenta of an event's leading charged leptons, whose distribution shapes display only a weak sensitivity to changes in √ s. This suggests the ability to build a search analysis resilient across varying collider energies. In Sec. 5, we turned to studying the behavior of hadronic observables over √ s = 14 − 100 TeV, for both the signal and background processes, and at various fixed-order and resummed precision. We focused exceptionally on hadronic activity in the context of a jet veto. When the veto's transverse momentum threshold (p Veto T ) is set to a static value of p Veto T = 30− 100 GeV for jet radii R = 0.1− 1.0 at √ s = 14− 100 TeV, we find discouraging survival efficiencies for signal benchmarks. Intriguingly, when using a dynamic jet veto, and in particular setting p Veto T on an event-by-event basis to the leading charged lepton p T in an event, we find qualitatively opposite behavior. Survival efficiencies for signal benchmarks reach ε(p Veto   T   ) 90 − 95% and remain largely independent of m N (for m N 200 GeV), R, and √ s; QCD scale uncertainties reduce at NLO+PS(LL) and NLO+NNLL for signal and background processes; and QCD background rejection capabilities in fact improve for increasing √ s. The impact of such a veto scheme on VBF-like topologies, events with finalstate τ leptons that decay hadronically, top quark and EW backgrounds, and, significantly, "fake lepton" backgrounds were also addressed. A summary of this work for a broad audience was first reported in Ref. [69].
Based on the findings of the previous two chapters, a new methodology is proposed in Sec. 6 to search for heavy neutrinos N in the trilepton final state, i.e., pp → N ℓ + X → 3ℓ + MET + X. The search analysis is based on employing a dynamic jet veto in conjunction with exclusive S T , which effectively discriminates according to the relative amount of leptonic and hadronic activities in an event. For various active-sterile mixing flavor hypotheses, the proposed analysis is compared with a state-of-the-art, benchmark analysis inspired by the analogous LHC Run II search [74]. At 14 TeV, we find that the proposed analysis can improve sensitivity to heavy Dirac neutrinos by an order of magnitude (in (m N , |V | 2 )-space) for m N 150 GeV. This improvement is achievable at both L = 300 fb −1 and 3 ab −1 , with the largest improvement occurring for the largest m N , and is independent of flavor hypothesis. In absolute terms, the TeV. See Sec. 6.3.3 for details. In none of the above cases was the proposed analysis optimized using multi-variant or machine learning techniques. We expect further improvement with such tools.
Our overarching computational framework, which largely takes advantage of out-ofthe-box, public Monte Carlo tools is described in great detail in Sec. 3. This includes the publication of a Dirac neutrino-variant of the FeynRules-based HeavyN libraries [62], and is available from the FeynRules database at feynrules.irmp.ucl.ac.be/wiki/HeavyN.
The anomalous production of charged leptons without hard, central hadronic activity in pp collisions is a key prediction of low-scale Type I Seesaw models. The prediction, however, is not unique: The success of the dynamic jet veto relies on the color structure of the hard, signal process and the ability to measure (even by proxy) the mass scale on an event-by-event basis. Hence, we believe much of our findings are immediately applicable to other searches for colorless particles that can be produced and decayed through colorless force carriers. For example: the production of exotically charged scalars and fermions, respectively from the Types II and III Seesaws; the production of heavy charged scalars that decay to long-lived stable, neutral particles, as commonly found in loop-level Seesaw and dark matter models; pair production of slepton and electroweakino pairs in Supersymmetry; and pair production of N through Z ′ gauge bosons in U (1) B−L theories. Moreover, we do not believe that setting the dynamic veto threshold p Veto T to the leading charged lepton in an event is always most optimal. (Indeed, we found in some cases that the subleading charged lepton gave better performance.) The exploration of other leptonic observables that are sensitive to the hard scale, e.g., S T and MET, is left to future studies. Similarly, vetoing on alternative hadronic observables, just at H T or jet mass, are also possibilities where further investigations are encouraged.

Conclusions
Colliders offer a complementary probe of the origin of active neutrino masses, their mixing, as well as their potential Majorana nature. If EW-or TeV-scale heavy neutrinos play a role in the generation of neutrino masses, and if such particles couple appreciably to the SM sector, then the LHC and its successors provide fertile ground for their discovery and, potentially, properties determination. In this study, we have systematically assessed the discovery potential of heavy neutrinos at pp colliders with √ s = 14 − 100 TeV.
This has been motivated by the unequivocally improved understanding of neutrino and jet physics, the sophistication of state-of-the-art Monte Carlo tools, and the mandates of the several HEP strategy updates ongoing at the time of this work. Relative to ongoing search-strategies [74], we find that the LHC's sensitivity can be improved by an order of magnitude in both the immediate term and over its lifetime. The increase in sensitivity can be attributed to a newly proposed search strategy (see Sec. 6) for heavy neutrinos decaying into a purely leptonic final state that, intuitively, relies on the ability to discriminate global leptonic activity from global hadronic activity. Crucially, this employs an unusual jet veto scheme, one where the p T threshold is set on an event-by-event basis to the p T of the leading charged lepton in the event. The functionality of this dynamic jet veto is described in technical detail in Sec. 5; a more broadly accessible summary is presented Ref. [69]. In addition, the proposed methodology exhibits by design an unusual resilience against variable collider energy, and therefore can serve as a baseline for future collider searches for multi-lepton searches of heavy neutrinos. The results of this study are encouraging and we look forward to the prospect of at last unveiling the mystery underlying tiny neutrino masses.