Majorana Neutrinos in Same-Sign $W^\pm W^\pm$ Scattering at the LHC: Breaking the TeV Barrier

We revisit the sensitivity to non-resonant, heavy Majorana neutrinos $N$ in same-sign $W^\pm W^\pm$ scattering at the $\sqrt{s}=13$ TeV LHC and its high-luminosity upgrade. As a benchmark scenario, we work in the context of the Phenomenological Type I Seesaw model, relying on a simulation up to next-to-leading order in QCD with parton shower matching. After extensively studying the phenomenology of the $pp\to\mu^\pm\mu^\pm j j$ process at the amplitude and differential levels, we design a simple collider analysis with remarkable signal-background separation power. At 95\% confidence level, we find that the squared muon-heavy neutrino mixing element $\vert V_{\mu N} \vert^{2}$ can be probed down to $0.1-0.5 ~ (0.05-0.1)$ for $m_N = 1-10~{\rm TeV}$ with $\mathcal{L}=300$ fb$^{-1}~(3$ ab$^{-1})$. For heavier masses of $m_N = 20~{\rm TeV}$, we report sensitivity for $\vert V_{\mu N} \vert^{2}\gtrsim 0.8~(0.4)$. The $W^\pm W^\pm$ scattering channel can greatly extend the mass range covered by current LHC searches for heavy Majorana neutrinos and particularly adds invaluable sensitivity above a few hundred GeV. We comment on areas where the analysis can be improved as well as on the applicability to other tests of neutrino mass models.


I. INTRODUCTION
Following the discovery of neutrino oscillations [1,2], uncovering the origin of neutrinos' tiny masses and their large mixing angles are among the most pressing questions in particle physics today [3,4]. To address these mysteries, neutrino mass models, collectively known as Seesaw models, do so by hypothesizing a variety of states that couple to the Standard Model's (SM) lepton and Higgs sectors [5]. These states include new charged or gauge-singlet (or sterile) fermions, scalars with exotic gauge quantum numbers, as well as gauge bosons of new symmetries, and mediate the non-conservation of lepton number (LN) and/or charged lepton flavor number over a range of mass scales and coupling strengths. For reviews of Seesaw models and their tests, see Refs. [6][7][8][9]. Despite these viable solutions, there remains a lack of clear guidance from both experiment and theory as to what is realized by nature. It is therefore necessary to broadly approach this challenge in complementary aspects.
To this extent, tests of neutrino mass models at the Large Hadron Collider (LHC) are supported by a number of signatures, including searches for dijet resonances [10,11], many-lepton final states [12][13][14][15][16][17] and LN-violating lepton pairs [13,[18][19][20], but rely mostly on mechanisms mediated by quark-antiquark (qq) annihilation [21]. However, due to its high center-of-mass energy * fuks@lpthe.jussieu.fr † jonas.neundorf@desy.de ‡ krisztian.peters@desy.de § richard.ruiz@uclouvain.be ¶ matthias.saimpert@cern.ch ( √ s), the LHC is also effectively an electroweak (EW) boson collider [22][23][24]. This in turn opens a multitude of complementary channels. For example, in the context of the Phenomenological Type I Seesaw model [6,25], the W γ fusion channel [26][27][28][29] has already helped direct searches for heavy neutrinos N with masses beyond a few hundred GeV improve sensitivity to active-sterile neutrino mixing matrix elements V N . In fact, with an integrated luminosity of L ≈ 36 fb −1 of proton-proton collisions at √ s = 13 TeV, |V N | 2 O(0.01 − 1) are excluded for lepton flavors ∈ {e, µ} and sterile neutrino masses in the range m N = 100 GeV − 1 TeV [13,18]. With the full LHC data set, this degree of sensitivity is anticipated to reach masses up to m N = 3 − 4 TeV [30]. Motivated by the recent experimental observations of EW vector boson fusion / scattering (VBF) at the LHC [31][32][33][34][35], we revisit the sensitivity of same-sign W ± W ± scattering to TeV-scale Majorana neutrinos at √ s = 13 TeV. As a benchmark scenario, we work in the framework of the Phenomenological Type I Seesaw and focus on the production of same-sign muon pairs (µ ± µ ± ) without substantial transverse momentum imbalance via the space-like exchange of an N in W ± W ± scattering [36], As shown in figure 1, this is essentially a realization of neutrinoless ββ (0νββ) decay at large momentum transfer when mediated at dimension d = 7 [37][38][39]. While past studies have investigated the importance of this channel [6,26,36,40,41], most works were restricted to sub-TeV m N , and therefore subject to signal processes with more dominant cross sections. As the heavy neutrino exchange in equation (1.1) is always non-resonant, the channel is complementary to other processes, such as the qq → N annihilation and W γ → N fusion mechanisms, which become kinematically inaccessible for sterile neutrinos that are too heavy. For similar reasons, the channel is robust against the impact of long N lifetimes.
To conduct this study, we employ a state-of-the-art simulation tool chain that allows us to model SM backgrounds and, for the first time, the W ± W ± → ± i ± j signal process at next-to-leading order (NLO) in QCD with parton shower (PS) matching. We report remarkable signal-background separation power and attribute this to the signal process exhibiting both VBF and LNviolating topologies. For the pp → µ ± µ ± jj collider signature with forward jet-tagging and simple selection cuts, we find that |V µN | 2 0.1−0.5 (0.05−0.1) can be probed at 95% confidence level (CL) for m N = 1 − 10 TeV with L = 300 fb −1 (3 ab −1 ). For masses of m N = 20 TeV we find sensitivity for |V µN | 2 down to 0.8 (0.4). The remainder of this work is organized in the following manner: First, we describe in section II our theoretical framework and give an overview of current experimental constraints. Next, we summarize our computational setup (section III) and simulation prescriptions (section IV). In section V we explore extensively the phenomenology of the W ± W ± → ± i ± j signal process at the amplitude and differential levels. This is then used in section VI to design our collider analysis and estimate the LHC's sensitivity to LN violation (LNV) in W ± W ± scattering. We provide an outlook in section VII on areas where our analysis can be improved as well as its applicability to other tests of neutrino mass models at colliders. Finally, we conclude in section VIII. Where relevant, technical derivations and details on software modifications are reported in appendices A, B, and C.

II. THEORETICAL FRAMEWORK
We describe in this section the theoretical framework in which we work (section II.1) and summarize current experimental constraints on the model (section II.2).

II.1. The Phenomenological Type I Seesaw Model
To study the sensitivity of the LHC to the LN-violating W ± W ± → ± i ± j process, we work in the context of the Phenomenological Type I Seesaw model [6,25]. Like the eponymous mechanism [43][44][45][46][47][48][49] or its low-scale variants [50][51][52][53][54][55], the model hypothesizes the existence of Majorana neutrinos (N k ) that couple to SM particles through mass mixing with light neutrinos (ν k ). More precisely, the renormalizable Lagrangian of the theory, extends the SM Lagrangian (L SM ) by kinetic and Majorana mass terms (L Kin. ) for n R ≥ 2 right-handed (RH) neutrinos (ν i R ), as well as by Yukawa couplings (L Y ) between the SM Higgs field, the SM lepton doublets (L j ) and the gauge-singlet fermions ν i R . After EW symmetry breaking, flavor eigenstates of active, left-handed (LH) neutrinos (ν L ) can be generically [6] decomposed into light (ν k ) and heavy (N k ) mass eigenstates with mass eigenvalues m ν k and m N k : Here, the complex-valued matrix elements U k (V k ) parameterize the mixing between the LH interaction state ν L and the light (heavy) mass eigenstate ν k (N k ). Formally, the matrix elements satisfy the relationship U U † + V V † = I. For clarity and without loss of generality, we order states such that m N k < m N k +1 . To leading order in active-sterile mixing |V k |, equation (2.2) gives rise to the following effective, charged current component of the theory's interaction Lagrangian, In the above expression, g W ≈ 0.65 is the SM weak gauge coupling constant and P L/R = (1/2)(1∓γ 5 ) are the standard chiral projection operators in four-component notation. Similar terms can be derived for the neutral current (Z) and Higgs interactions for both heavy Dirac and Majorana neutrinos [6]. For simplicity we solely investigate the phenomenology of the lightest heavy neutrino mass eigenstate, which we relabel as N ≡ N k =4 (so that V N ≡ V k =4 ), and decouple all other heavy eigenstates.
In short-distance scattering and decay processes involving only a single heavy neutrino, the mixing factors that appear in equation (2.3) act as effective couplings and factor out of amplitudes. For resonant production of a heavy neutrino, this allows one to define a "bare" cross section σ 0 such that [56] σ(pp → N ± + X) ≡ |V N | 2 × σ 0 (pp → N ± + X). (2.4) When the W ± W ± → ± i ± j process involves only a single t-channel exchange of a heavy neutrino, such a factorization is also possible. More specifically, one can define a "bare" t-channel cross section given by For our purposes, it suffices to add that these expressions hold at NLO in QCD [29,57]. Under certain assumptions, the above decomposition or a similar one can hold for processes involving multiple interfering heavy neutrino mass eigenstates. Such expressions, however, may not always be tractable due to large interference effects. For further discussions, see Refs. [58][59][60][61][62].

II.2. Model Constraints
In its most general construction, the free parameters of the Phenomenological Type I Seesaw model (beyond those of the SM) consist of the neutrino masses m ν k , m N k , and the neutrino mixing elements U k , V k . Imposing flavor symmetries on the lepton sector, however, can reduce the number of independent degrees of freedom, as discussed in Refs. [58-60, 63, 64] and references therein. It is also possible to constrain these parameters by tying the lightness of the ν k neutrinos to beyond the SM physics, such as to dark matter and the baryon asymmetry of the universe, as done for example in Refs. [65,66]. For our purposes, we take mass and mixing parameters to be phenomenologically independent. We do this to develop a collider analysis in a flavor-model-independent fashion and is motivated by the desire to broaden sensitivity to a range of ultraviolet completions.
Beyond theoretical considerations are the following recent experimental constraints on the model: • Direct constraints from 0νββ searches: After an exposure of 127.2 kg-yrs, the GERDA experiment reports at 90% CL the following lower limit on the 0νββ decay half-life in 76 Ge [67]: Assuming that the nuclear 0νββ process is only mediated by heavy neutrinos, this translates into an upper limit on their masses and mixing of where the variation stems from the uncertainties in the nuclear matrix element. For further details on the derivation of this constraint, see appendix A.
• Direct constraints from collider searches: At √ s = 13 TeV and with L ≈ 36 fb −1 of data, searches for the pp → i j k + E miss T signature with ∈ {e, µ}, by the CMS experiment constrain active-sterile neutrino mixing at 95% CL to be [13], Constraints from the ATLAS experiment with the same integrated luminosity are comparable for m N < m W but weaker for m N > m W due to the absence of the W γ channel in their signal modeling [15]. Searches for the pp → ± i ± j + nj signature by CMS yield only slightly more stringent constraints due to a larger signal-over-background ratio [18].
• Indirect constraints on V k : For n R = 3 sterile neutrinos with masses above the EW scale, a global study of precision EW data, searches for nonunitarity in quark mixing, and searches for lepton flavor violation and non-universality of weak decays constrain active-sterile neutrino mixing to be [68], at 95% CL. The parameter η is related to the heavy neutrino mixing matrix V k by 2|η | = For the scenario that we consider, i.e., = = µ with only one heavy neutrino species, this translates to an upper limit of |V µN | 2 < 4.41 · 10 −4 at 95% CL . (2.10) • Direct constraints on the absolute mass scale of light neutrinos: Attempts to measure the light neutrino mass scale directly from the kinematic end point in β decay with the KATRIN experiment [69] constrain the light neutrino masses to satisfy [70] m(ν e ) < 1.1 eV at 90% CL . (2.11) • Constraints on neutrino masses from cosmology: An analysis of neutrinos' impact on the cosmic microwave background, supernovae, large scale structure, and big bang nucleosynthesis constrains the sum of light neutrino masses to be [71] m ν 0.26 eV at 95% CL . (2.12) • Neutrino oscillation measurements of U k : In the absence of sterile neutrinos or additional new physics, the elements of the light neutrino mixing matrix U k have been fit to or constrained by long and short baseline neutrino oscillation data under the condition that U is unitary [72]. Relaxing this constraint, however, greatly weakens the goodness of fit, particularly in the τ flavor sector [73,74].

III. COMPUTATIONAL SETUP
We now summarize the computational setup used in this study. We start with section III.1 where we document our Monte Carlo (MC) simulation chain. In section III.2 we list the numerical values used for SM inputs, and similarly in section III.3 the numerical values used for non-SM inputs. A description on how we model the W ± W ± → ± i ± j signal process and leading background processes in pp collisions is deferred to section IV.

III.1. Monte Carlo Setup
To investigate same-sign W ± W ± scattering when mediated by a heavy Majorana neutrino at the √ s = 13 TeV LHC, we employ a state-of-the-art simulation tool chain. For hard, parton-level scattering processes, we use the general-purpose MC event generator Mad-Graph5 aMC@NLO (MG5aMC) (version 2.7.1.2) [75,76], which enables us to simulate tree-induced processes in the SM up to NLO in QCD [77][78][79][80]. Processes involving heavy neutrinos are simulated up to NLO in QCD by importing into MG5aMC the default variant of the HeavyN [29] FeynRules UFO libraries [81][82][83][84][85]. For select backgrounds we perform jet matching up to one additional parton (at the Born level) at NLO in QCD precision using the FxFx matching procedure [86], as further detailed in section IV. To simulate the decay of W bosons, we impose the spin-correlated narrow width approximation as implemented in MadSpin [87,88].
All signal and background events are passed through Pythia8 (version 243) [89] for QCD and QED parton showering up to leading logarithmic (LL) accuracy, hadronization, and multiparton interaction / underlying event modeling. Decays of heavy-flavored hadrons and τ leptons are handled internally by the PS. Following Refs. [90,91], we apply the improved color reconnection [92] and dipole recoil [93] models available in Pythia8. For simplicity, activity from additional proton-proton interactions occurring during the same bunch crossing, i.e., pileup, is assumed to be subtracted from experimental data by dedicated algorithms and, therefore, are not included in our simulations.
In our analysis, hadron-level events are passed to the fast detector simulator Delphes (version 3.4.2) [94]. There, particle-level clustering of hadrons is handled according to the anti-k T algorithm [95][96][97] as implemented in FastJet [98,99], with a radius parameter R = 0.4. To emulate experimental reconstruction with realistic detector resolution and particle identification, detector responses are tuned using the ATLAS configuration card available in the Delphes repository. Particle-level distributions at LO+PS and NLO+PS were checked using MadAnalysis5 (version 1.8) [100][101][102].

III.2. Standard Model Inputs
For SM inputs we work in the n f = 4 active quark flavor scheme with a Cabibbo-Kobayashi-Maskawa mixing matrix that is diagonal with unit entries. Unless specified, we assume the following mass and coupling values: For scattering computations at both LO and NLO in QCD we employ the NNPDF3.1 NLO+LUXqed parton distribution function (PDF) set with α s (M Z ) = 0.118 (lhaid=324900) [103][104][105]. For non-perturbative dynamics we tune the shower with the ATLAS A14 central tune (Tune:pp=21) as paired with the NNPDF2.3 LO+QED PDF set with α s (M Z ) = 0.130 (pdfcode=247000) [106]. PDF and α s scale evolution are handled by LHAPDF (version 6.2.3) [107]. PDF uncertainties are extracted using replicas PDFs [105,107]. For signal and background processes, we set the nominal (ζ = 1.0) collinear factorization (µ f ) and QCD renormalization (µ r ) scales to be half the sum over each visible, final-state particle's transverse energy: Here m f and p T,f stand for the mass and transverse momentum of the final-state particle f respectively. The shower factorization scale (µ s = ζ ×μ s ) is kept at its default value prescribed in Ref. [76]. To estimate the size of higher-order QCD corrections, we vary discretely and independently µ f , µ r , µ s , over the set ζ = {0.5, 1.0, 2.0} to obtain a nine-or 27-point uncertainty band.

III.3. Heavy Neutrino Inputs
For simulations involving the heavy neutrino N we use the default inputs of the Majorana neutrino variant of the HeavyN NLO UFO libraries [29]. SM particle masses are updated according to equation (3.1). Additional heavy mass eigenstates are decoupled by setting m N5 , m N6 = 10 10 GeV. To deactivate e and τ flavor mixing, we set Sensitivity to smaller values of |V µN | are obtained by a naïve rescaling of cross sections, which is permissible by equation (2.5). Notably, as N is never resonantly produced, its total width (Γ N ) and lifetime can be ignored.

IV. SIGNAL AND BACKGROUND MODELING IN MADGRAPH5 AMC@NLO
Generically speaking, collider processes that feature either LNV or VBF exhibit characteristic kinematical and topological properties that enable remarkable background rejection capabilities. Consequentially, background modeling for the W ± W ± → ± i ± j channel is hindered by high background rejection rates, and hence by poor MC efficiencies. In this context, we report the development of efficient MC modeling prescriptions that overcome such difficulties for our signal (section IV.1) and leading backgrounds (section IV.2) at NLO+PS within the MG5aMC simulation framework.

IV.1. Signal Modeling
To model the W ± W ± → ± i ± j process when mediated by the exchange of a t-channel Majorana neutrino N , as shown in figure 1, we consider the gauge-invariant set of W ± W ± scattering diagrams contributing to the following 2 → 4, Born-level process at O(α 4 ) Here q denotes any light quark or antiquark. We neglect interference with the s-channel process, To justify this, we require that the leading dijet system at the analysis level carries a large invariant mass, which suppresses the N ( * ) → W ( * ) → qq splitting-chain. For m N below the TeV scale, neglecting the qq annihilation mechanism is also justified by the narrow width approximation. Under this the leading contributions to the s-channel process are factorizable and non-interfering with t-channel diagrams; interfering contributions are O(Γ N /m N ) 1, and hence insignificant. To carry out this modeling in MG5aMC at NLO in QCD using the HeavyN libraries, we employ the syntax 1 import model SM_HeavyN_NLO define p = g u c d s u~c~d~sd efine j = p generate p p > mu+ mu+ j j QED=4 QCD=0 $$ w+ w-/ n2 n3 [QCD] add process p p > mu-mu-j j QED=4 QCD=0 $$ w+ w-/ n2 n3 [QCD] Formally, the Born-level matrix element for equation (4.1) is finite in the absence of phase space cuts. At O(α s ), however, an infrared-safe definition for external states is needed. We therefore require at the generatorlevel that QCD partons are sequentially clustered according to the anti-k T algorithm with R = 0.4 and that 1 For further details on syntax and usage, see Refs. [29,76].
the transverse momentum (p T ) and pseudorapidity (η) of these clusters satisfy the following demands: p j T > 20 GeV and |η j | < 5.5. (4.2) As a technical remark, we relax checks on infrared pole cancellation in MadFKS. Such checks are automatically raised in MG5aMC for VBF processes due to the possible omission of virtual diagrams that are mixed NLO QCD-EW corrections and not pure QCD contributions at O(α s ). In our case such diagrams do not exist and therefore bypassing the check only impacts the MC efficiency. To do this we set #IRPoleCheckThreshold=-1.0d0 in the file Cards/FKS params.dat.

IV.2. Background Modeling
Due to the presence of forward, high-p T jets, a highenergy, same-sign lepton pair, and an absence of (light) neutrinos in our signal process, several backgrounds processes that are traditionally present in collider searches for LNV can be readily suppressed through simple kinematic requirements. Among many examples are: demanding that the leading dijet system carries a large invariant mass, stringent p T cuts on the same-sign leptons, and vetoing events with three or more charged leptons. As a result, background categories such as associated top quark production, ttB with B ∈ {W ± /γ * /Z, h}; single top quark channels, tB; and triboson production W W V can be neglected for the purposes of our study.
The leading background channels that remain after such baseline selection criteria include: the mixed EW-QCD channel W ± W ± jj (section IV.2.1), the pure EW channel W ± W ± jj (section IV.2.2), and the mixed EW-QCD diboson+jets process W ± V + nj with V ∈ {γ * /Z} (section IV.2.3). To model these processes in MG5aMC, we employ the prescriptions described below.
IV.2.1. QCD Production of Same-Sign W ± W ± jj Due to its similar topology and kinematic scales, the mixed EW-QCD production (henceforth labeled QCD production) of pp → W ± W ± jj is a prominent background for the pp → ± ± jj signal process. However, a defining characteristic of this mode, which at the Born level occurs at O(α 2 α 2 s ), is the t-channel exchange of a gluon. This indicates an intermediate flow of color. As such, the presence of a third, high-p T jet is significantly more likely in this process than in the signal process. Such radiation can induce recoils in the (W ± W ± )system, siphon energy from the two forward jets, or give rise to excess central hadronic activity [108][109][110][111][112]. One is thus motivated to consider the process at NLO in QCD.
We do this in MG5aMC by using the syntax:  add process p p > w-w-j j QED=2 QCD=2 [QCD] As in the signal process, we impose the criteria of equation (4.2) on outgoing QCD partons. We report in the first line of table I the corresponding generator-level cross section at √ s = 13 TeV, which reaches σ ∼ 385 fb, along with residual µ f , µ r scale and PDF uncertainties. Before parton showering, resonant W bosons are decayed to muons (see section III.1 for related details). For simplicity, we neglect contributions from leptonic tau decays. We do so because the presence of additional light neutrinos results in events with characteristically softer muons and larger momentum imbalances. Such features can be tamed by tuning the pre-selection and signal region cuts used in section VI.

IV.2.2. EW Production of Same-Sign
Like the previous case, pure EW production of pp → W ± W ± jj, which at the Born level occurs at O(α 4 ), is a leading background to the pp → ± ± jj signal process. Similarities to the signal process's topology, kinematic characteristics, and, importantly, color flow are identifiable in the W ± W ± → W ± W ± , VBF sub-process. While featuring a large cross section at the inclusive level, other sub-processes, such as resonant triboson production, are essentially removed through basic selection criteria.
To efficiently model the EW same-sign W ± W ± channel at NLO in QCD, we make a variant of the so-called "vector boson fusion approximation" and neglect resonant triboson production in a gauge-invariant manner. We do this under the presumption that relevant analysislevel selection cuts are applied. Other non-resonant, interfering diagrams are kept. The syntax we employ is generate p p > w+ w+ j j $$ a z w+ w-QCD=0 QED=4 [QCD] add process p p > w-w-j j $$ a z w+ w-QCD=0 QED=4 [QCD] We apply the same loose, generator-level cuts as listed in equation (4.2) and report good numerical stability. In the second line of table I, we report the generator-level cross section at √ s = 13 TeV, which reaches σ ∼ 254 fb, and associated uncertainties. Resonant W bosons are decayed to muons before parton showering.

IV.2.3. Inclusive Diboson Spectrum
Due to its picobarn-scale rate, the inclusive pp → 3 ν + X spectrum (which we imprecisely label "diboson spectrum" even though the process includes interference with all non-resonant diagrams) contributes to the same-sign dilepton signature pp → ± ± jj through pathological configurations of the final-state kinematics. Such configurations include, for example, when two or more initial-state QCD emissions both possess large p T but the odd-sign charged lepton is not successfully identified because it is too soft in p T or too forward in |η|.
While a bulk of these phase space configurations are captured by the fixed-order matrix element for the pp → 3 νjj process, which at the Born level occurs at O(α 4 α 2 s ), convergence in perturbative QCD requires that outgoing QCD partons are hard (high p T ) and central (low |η|). In light of the anticipated scales of the (3 ν)-system, the component of phase space where one QCD parton is central and one is forward is better described by the fixedorder matrix element for the pp → 3 νj sub-process, with the second j being populated by the PS. Likewise, for two forward emissions, the phase space is better described best by the pp → 3 ν sub-process with two PS emissions.
To model these complications we extend the bruteforce prescription of Ref. [30]. This entails starting with the 2 → 4 process, pp → 3 ν at NLO in QCD, with all interfering, non-resonant diagrams and without invoking the narrow width approximation for intermediate W, Z bosons. Instead, loose, generator-level cuts on leptons are imposed to regulate s-and t-channel divergences, thereby keeping matrix elements finite and perturbative. Specifically, we require that charged leptons satisfy m os > 8 GeV and |η | < 4.0, where m os is the invariant mass of any opposite-sign, charged lepton pair, independent of flavor. While smaller m os thresholds can still regulate γ * → splittings, we refrain from doing so to avoid contributions from vector meson resonances. Such states are not modeled with the perturbative event generator MG5aMC.
To account for additional central jet multiplicities, we match the inclusive pp → 3 ν spectrum at NLO+PS up to its first jet multiplicity (relative to the Born level) at NLO in QCD using the FxFx prescription [86]. The relevant MG5aMC syntax in this case is define ell = e+ mu+ ta+ e-mu-tadefine vv = ve vm vt ve~vm~vtg enerate p p > ell ell ell vv QED=4 QCD=0 [QCD] @0 add process p p > ell ell ell vv j QED=4 QCD=1 [QCD] @1 Explicitly, we use the following jet-matching inputs: where Q FxFx cut is the FxFx matching scale. With this setup, hadronic observables are accurate to at least [fb] Upper: As a function of heavy neutrino mass mN [GeV], the bare cross section σ/|V N V N | 2 [fb] for the W ± W ± signal process at NLO in QCD (purple band), as well as the bare cross sections σ/|V N | 2 for the CCDY (black band) and W γ fusion (green band) processes at NLO in QCD. Band thickness corresponds to the residual scale uncertainty. Lower: The QCD K-factor for each channel.
LO+PS(LL), with the two leading jets being defined at all momenta and rapidities. We report in the third line of table I the generator-level cross section of the diboson spectrum at NLO in QCD with FxFx matching to the first jet multiplicity (FxFx1j) and its associated uncertainties. At √ s = 13 TeV the rate is about σ ∼ 2.5 pb, which is just slightly larger than the rate at NLO, which we compute to be σ ∼ 2.3 pb.
With this setup, we capture configurations where the odd-sign charged lepton in the 3 ν2j +X final state is too forward or too soft to be identified as an analysis-quality charged lepton. A disadvantage of this setup, however, is the limited MC statistics when the two same-sign charged leptons carry p T 100 − 150 GeV but the odd-sign lepton is much softer. To enrich MC statistics for this region of phase space, we introduce tailored generator-level cuts into the MG5aMC phase space integration routines. Enriched samples are combined with the baseline FxFx1j sample. Overlap is removed through cuts on p 2 T . For technical details of this modeling, see appendix B.
In this section we investigate the phenomenology of the W ± W ± → ± i ± j process when mediated by a heavy Majorana neutrino at the LHC. To do this, we examine the integrated (section V.1) and differential (section V.2) cross sections of the W ± W ± channel, and place special focus on the low-(section V.1.1) and high-mass (section V.1.2) limits of the intermediate neutrino, on the impact of QCD corrections (section V.1.3), and on potential violations of partial-wave unitarity (section V.1.4).

V.1. Total Production Rate
As a first step, we present in the upper panel of figure 2 and as a function of heavy neutrino mass m N , the total cross section for the full 2 → 4, hadron-level process pp → ± ± jj + X.
For the mass range m N = 40 GeV − 20 TeV, we report that bare cross sections at NLO in QCD, QCD K-factors, as well as scale and PDF uncertainties roughly span A summary of bare cross sections and uncertainties for the W ± W ± → ± ± process at representative heavy neutrino masses is listed in table II. To compare to other heavy neutrino processes, 2 we also present in figure 2 the bare cross sections at NLO in QCD, as defined in equation (2.4), the associated scale uncertainties, and QCD K-factor for the 2 → 2, charged current Drell-Yan (CCDY) process (black band), and the 2 → 3, W γ fusion process (green band), We find several notable observations: First is that, quantitatively, the bare, same-sign W W cross section is about 4 − 6 (1 − 3) orders of magnitude smaller than the CCDY (W γ) process for m N ∼ 50 − 100 GeV. This is much smaller than the 4 (2) orders of magnitude that one expects from naïve power counting. Second is that while the bare rates of resonant channels fall precipitously for increasing m N , which is due to suppression in both the matrix element and available phase space, the W ± W ± rate moderately increases before slowly decreasing. For an active-sterile mixing of |V N | 2 = 1, this leads to the W ± W ± rate surpassing the W γ rate at m N ∼ 500 − 600 GeV and the CCDY rate at m N ∼ 700 GeV. Due to the different sensitivities of the three channels to active-sterile mixing, the crossover occurs at higher neutrino masses for smaller values of |V N | 2 . For example: While mixing can alter the precise values of these crossovers, the qualitative picture does not change. For instance: independent of |V N | 2 , the W ± W ± → ± ± cross section exhibits a qualitatively different dependence on m N than in the CCDY and W γ channels. This leads to the W ± W ± rate at √ s = 13 TeV to be the same at both m N ∼ 40 GeV and m N ∼ 2.5 TeV. Moreover, unlike resonant production of heavy neutrinos via qq annihilation or W γ fusion, heavy neutrinos in t-channel processes like W ± W ± → ± i ± j are non-resonant. So while there is a kinematic suppression in the W ± W ± matrix element at very large m N , there is no corresponding phase space suppression. This manifests in the cross section as a milder dependence on increasing sterile neutrino masses. Interestingly, as heavy neutrinos are never on-shell in t-channel exchanges they can never manifest as a long-lived particle. Therefore, search complications associated with displaced vertices are not present.
In comparison to past work, this is the first evaluation of the full 2 → 4, same-sign W ± W ± scattering pro-cess in pp collisions at NLO in QCD. At LO, the literature [6,36,39,40] is admittedly in disagreement with itself. Qualitatively, the dependence on collider energy and heavy neutrino masses in all these works are consistent. Quantitatively, large differences exist. In some cases, differences can be traced to the omission of numerical pre-factors in analytic and/or numerical results, theoretical uncertainties associated with the effective W approximation [22][23][24], and uncertainties in PDF sets. In other cases, the lack of documented inputs and possible phase space cuts hinder precise comparisons. Support for our numerical results include agreement with analytical expressions for helicity amplitudes. For further details, we refer to sections V.1.1 and V.1.2, and appendix C.
To further understand the dependence on m N in the W ± W ± channel, we consider for illustration purposes the matrix element for the 2 → 2, W ± W ± → ± ± subprocess. For the momentum and helicity assignments Here ε are the usual helicity polarization vectors for massive gauge bosons in the unitary gauge, the (t ↔ u) term accounts for final-state lepton exchange, and following the Feynman rules of Refs. [113,114], the LN-violating current (T ) in the HELAS convention [115] is 3 [36] T In the above we assume a clockwise fermion flow of leptons [113,114], p N = (p W 1 − p 1 ) is the momentum of the internal sterile neutrino, and u(p, λ) and v(p, λ) are the standard helicity spinors for massless, spin-1/2 fermions. Crucially, differences in Feynman rules for LNviolating fermion currents relative to the standard rules for LN-conserving ones give rise to an effective parity inversion in the W + 1 − + 1 − N vertex and spinor for + 1 [116,117]. This implies [118,119] that the successive gauge interactions involving massless particles in T are helicity inverting and not helicity preserving as one usually finds in SM gauge interactions involving massless, external particles. Subsequently, projection operators select for the heavy neutrino's RH helicity state, and hence the factor of m N × I 4 in T . This is in contrast to LN-conserving currents, such as in the process W + W − → + − , where projection operators select for the LH helicity state, and hence the p N term in T .
We report that exact, analytic evaluation of equations (5.10) and (5.12) yields somewhat bulky expressions without obvious insights. This is despite being a 2 → 2 process and can be tied to the added algebraic complication of the incoming W bosons being massive and carrying a longitudinal polarization. Instead, we focus on the low-(section V.1.1) and high-mass (section V. We consider first the limit where masses of both W and N are small compared to the W ± W ± scattering scale, i.e., when m W , m N M W W . In this limit, the LNviolating tensor current in equation (5.12) scales as where the M W W factor in the numerator originates from the two lepton spinors, u(p ), v(p ) ∼ √ M W W . In this same limit, the scattering of longitudinally polarized W bosons is enhanced over the scattering of the transverse polarizations. This enhancement can be seen in the polarization vectors themselves, which scale as This shows that in the high-energy limit the W ± W ± → ± ± process is driven by W ± 0 W ± 0 scattering and that the corresponding matrix element scales as Remarkably, after squaring M, the quadratic dependence on M W W is canceled by the flux factor in the definition of the parton-level cross sectionσ. This renders the rate independent of M W W but quadratic in m N , Thus, we can attribute the growth in the same-sign W W scattering rate seen in figure 2 for sub-TeV heavy neutrinos to the cancellation of momentum scales in highenergy W ± 0 W ± 0 scattering in tandem with helicity inversion in the LN-violating lepton current.
After a more careful computation (see appendix C), the W + W + → + 1 + 2 cross section for n R heavy neutrinos iŝ

. High-Mass Limit
We consider now the kinematic limit where the W boson's mass and all momentum-transfer scales are small compared to the sterile neutrino's mass, i.e., the decoupling limit [120] where m 2 W , M 2 W W , |t|, |u| m 2 N . In this limit, N -exchanges can be treated as contact interactions, and the pole structure of its propagator can be systematically expanded. Doing this causes the LN-violating tensor current of equation (5.12) to scale as An analogous expression holds for the u-channel.
A consequence of this expansion is that the angular dependence that encapsulates forward-and backwardscattering enhancements in gauge interactions becomes a sub-leading contribution in the propagator. This implies that for most W ± W ± polarization combinations the forward (t-channel) and backward (u-channel) helicity amplitudes are indistinguishable, save for a relative minus sign that triggers an exact destructive interference. As a result, the only non-vanishing amplitudes are those where the incoming W ± W ± states carry the same helicity.
Noting once more the enhancement of longitudinallongitudinal scattering over other W ± W ± helicity configurations, one finds that the leading contribution to the matrix element for the W ± 0 W ± 0 → ± ± process scales as After squaring, the dependence on M W W is partially compensated by the flux factor in the parton-level cross section. The result is a total rate that scales aŝ Immediately, we see that the m −2 N factor originating from the heavy neutrino propagator is never fully offset by the LN-violating current or other factors at the cross section level. Ultimately, this is responsible for the drop in cross section that occurs in figure 2 for increasing m N .
After a more careful computation (see appendix C), one finds that the parton-level, W + W + → + 1 + 2 cross section for the exchange of n R heavy neutrinos is: We stress that the transition rate's dependence on masses and mixing elements of heavy neutrinos mirrors the scaling behavior found in nuclear 0νββ decay rates [6,36].
Notably, the coherent summation over V permits complex phases to trigger potentially large cancellations in analogy to the "funnel behavior" in 0νββ decay.
Returning to figure 2, we recall that past investigations into the W ± W ± → ± i ± j process historically [36,40,41] employed the effective W approximation [22][23][24]. In this approximation, W bosons are treated as constituents of the proton and the 2 → 2, W ± W ± → ± i ± j scattering rate is convolved with W boson PDFs of the proton. Only later in Ref. [6] was the full 2 → 4 process with forward jets evaluated. In all these cases, however, only LO estimates of cross sections were calculated. Therefore, we are motivated to comment on the size of NLO in QCD corrections and residual uncertainties in the full 2 → 4 process, particularly in relation to those of the CCDY and W γ modes, which were first reported in Refs. [29,121].
In the lower panel of figure 2 we show as a function of heavy neutrino mass the QCD K-factors, as defined in equation (5.2), for the same-sign W ± W ± , CCDY, and W γ processes. Band thicknesses correspond to the residual µ r , µ f dependence at NLO. Over the mass range m N ∼ 40 GeV − 20 TeV, we find that QCD corrections gradually and uniformly increase the total W ± W ± rate by +5% to +35%. As in deeply inelastic scattering, one-loop QCD corrections to space-like EW emissions in VBF do not appreciably alter cross section normalizations [122]. Consequentially, we attribute the increase in cross section at NLO to real, initial-state radiation. This purported reliance on an O(α s (µ r )) splitting is supported by the increased scale dependence at larger m N and the lack of a scale dependence in the 2 → 4 process. In comparison to the other channels, the K-factor for W ± W ± sits just below (above) the CCDY (W γ) curve for m N ∼ 750 GeV and overtakes both at m N 1 TeV.
For corrections beyond NLO in QCD, one can consider two complementary directions. The first pertains to higher-order QCD corrections while the second pertains to EW corrections. Based on results for VBF production of the SM Higgs boson [123][124][125], we anticipate that improvements at O(α 2 s ) and O(α 3 s ) have only a modest impact on total cross sections and distributions. At NLO in EW, corrections to the LO rate typically scale as As we show in section V.2, scales for M W W at the LHC range about M W W ∼ 300 GeV − 600 GeV for a large array of heavy neutrino masses. This translates to a modest uncertainty of δσ NLO−EW /σ LO ∼ 9% to 13%. We anticipate that in both cases the impact of missing higher-order corrections is negligible for discovery purposes.
As a brief remark, we interestingly note that the scattering amplitude for the W ± W ± → ± ± process exhibits poor high-energy behavior when initiated by a pair of longitudinally polarized W bosons. As evident in equations (5.19) and (5.25), and more precisely in appendix C, matrix elements for both the low-and high-mass limits scale with some positive power of the (W W )-scattering energy, M W W . This is distinct from the CCDY channel where no such scaling behavior is present. For fixed heavy neutrino masses and mixing, such a dependence on M W W implies [126][127][128][129][130][131] that the matrix elements violate partial-wave unitarity above some scattering energy threshold E U , unless additional physics cancels this dependence. As both amplitudes depend on the mass of an internal Majorana neutrino, it is possible that the unitarity violation is actually tied to the explicit breaking of LN symmetry in the Type I Seesaw model. If so, then it can potentially be resolved through the spontaneous breaking of LN symmetry via a Higgs-like mechanism.
While a systematic study of partial-wave unitarity in the Phenomenological Type I Seesaw model lies beyond the scope of this work, we can nevertheless provide a qualitative outlook. Following Ref. [120], the J = 0 partial-wave amplitude of the W ± 0 W ± 0 → ± ± process is related to its matrix element M by the relationship Here θ 1 is the polar angle of 1 in the frame of the (W W )system. For the low-and high-mass limits, the partialwave amplitudes are given to lowest order by low-mass : a 0 ≈(2 − δ 1 2 ) Requiring |a J | < 1 implies that the W ± 0 W ± 0 channel saturates unitarity at M W W = E U , with low-mass : In the low-mass limit and assuming a single heavy neutrino of mass m N = 1 TeV with an active-sterile mixing of |V N | 2 = 10 −1 (10 −2 ) [10 −3 ], partial-wave unitarity saturates at about E U ∼ 7.6 TeV (76 TeV) [760 TeV]. Working instead in the high-mass limit, the mass and mixing upper bound from equation (2.7) as derived from 0νββ decay searches implies a lower bound on the saturation scale of E U ∼ 72 TeV − 87 TeV. This suggests that if a discovery of 0νββ decay is made by currentgeneration experiments, and if the decay is mediated by a heavy Majorana neutrino, then a future pp collider with √ s = 100 TeV may be able to probe this high-energy behavior. However, we caution that in both cases the saturation scale is acutely sensitive to the values of activesterile mixing elements and heavy neutrino masses.
We now turn to exploring the kinematic properties of the W ± W ± → ± i ± j signal process at √ s = 13 TeV. As NLO-vs-LO comparisons of VBF kinematics are extensively documented, we restrict ourselves to NLO+PS(LL) distributions where available and neglect comparisons to properties at LO+PS(LL). For concreteness, we fix i = j = µ and set simulation inputs as prescribed in section III. For each of the following observable we assume the representative benchmark masses m N = 750 GeV (darkest), 1.5 TeV (dark), and 5 TeV (light). Events are normalized to L = 300 fb −1 . Also shown for each distribution is the residual µ r , µ f , µ s uncertainty (band thickness) as obtained from a 27-point variation envelope. Throughout this section we work with particle-level objects to emulate detector thresholds (but not detector resolution) and ensure the infrared safety of observable definitions. This means that after parton showering we impose anti-k T (R = 0.4) clustering on all hadronic activity. We also require that electron, muon, hadronic tau,  -102]. That is to say, we consider an ideal setup in which identification efficiencies are set to unity and mistagging rates are set to zero. No kinematic smearing is applied. Given these stipulations, we define our signal as the pp → µ ± µ ± jj + X (5.36) process, where X denotes the possibility of additional hadronic or photonic activity. More precisely, we require events to possess exactly two same-sign µ candidates and at least two j candidates. Events containing any number of e or τ h candidates are rejected. This setup implies that we remain inclusive with respect to soft and forward objects that fail candidacy requirements. For clarity, objects are ranked by their p T , with p k T > p k+1 T . We start with figure 3 where we plot the (a,b) p T and (c,d) η distributions of the (a,c) leading (µ 1 ) and (b,d) sub-leading (µ 2 ) muon in our signal process at NLO+PS. We foremost note the lack of any resonant structure in both p µ T distributions. This follows from the absence of schannel resonances in the W ± W ± → µ ± µ ± sub-process. Instead, we find a p T behavior reminiscent of open particle production that plateaus for a few hundred GeV and then falls due to kinematic suppression. There is a steeper falloff for smaller m N . The p T spectra indicate that the (µ ± µ ± )-system, or equivalently the (W ± W ± )system, possesses a large invariant mass that reaches several hundred GeV. We observe that both muons tend toward smaller values of |η|, independent of the heavy neutrino mass, indicating an absence of forward scattering.
Moving onto figure 4, we present the same information as in figure 3 but for the (a,c) leading (j 1 ) and (b,d) subleading (j 2 ) jets. In both p T spectra we observe peaks at p j T ∼ m W /2, which is characteristic of the VBF process and is due to the recoil against the t-channel emission of W bosons. Consistently, we observe in the η distributions that the two jets are forward, with maxima in the forward direction near |η j | ∼ 3 and a suppression of central activity at |η j | ∼ 0. We find that the shapes of all observables in figure 4 are insensitive to the values of m N under consideration. This follows from the fact that quarks in the 2 → 4 process do not directly couple to the LN-violating current, and therefore act like spectators. As the pp → µ ± µ ± jj + X process is simulated at NLO+PS, one has access to the pp → µ ± µ ± jjj + X process at LO+PS. We are thus able to explore the QCD radiation pattern in the VBF process. In figure 5 we show the (a) p T and (b) η distribution of the trailing jet j 3 at LO+PS. For all considered values of m N , we observe an inclination toward lower p T , with most of the phase space sitting between the threshold at p j3 T ∼ p T = 25 GeV and p j3 T ∼ 40 GeV. This is just below the characteristic p T of the two leading jets. We note a strong suppression of central (|η j3 | 2) tertiary jets. This does not mean an absence of QCD radiation for |η j3 | 2, only that it is soft. Most of the activity resides in the forward direction, peaking at |η| ∼ 3−4, again independent of m N . As this is well in the vicinity of the leading jets it is likely that the q → qg and g → qq splittings responsible for j 3 involve smaller momentum transfers, which results in shallower opening angles between j 3 and its companion.
To further investigate the dynamics of the W ± W ± → µ ± µ ± sub-process, we consider in figure 6 observables that are built from the momenta of two or more particles. We start with figure 6(a), where we plot the azimuthal separation of the same-sign µ ± µ ± pair, defined as We find that the leptons exhibit a strong back-to-back trajectory with the separation peaking (curtailing) at ∆ϕ(µ 1 , µ 2 ) ≈ π (0). This is despite being a 4-body final state at LO, which would suggest a sizable recoil against the (jj + X)-system. For increasing heavy neutrino masses we observe a higher tendency for back-toback trajectories. The marginal-to-moderate recoil that is found suggests that modeling the 2 → 4 signal process as a 2 → 2 process within the effective W approximation as done in Refs. [36,40,41] is a fair approximation. In figure 6(b) we focus on the distribution of the missing transverse energy E miss T , defined per event as the magnitude of the two-momentum recoil against all vis- We find that the distribution strongly peaks at E miss T 10 GeV, in line with expectations of a 2 → 4 process without outgoing light neutrinos. As we are working without any detector resolution effects, the nonzero E miss T originates from the weak decays to light neutrinos of mesons generated in the parton shower. Aside from differences in the rate normalization, we observe no substantial dependence of E miss T on the heavy neutrino mass. In figure 6(c) we show the invariant mass distribution of the two highest p T (leading) jets, given by M (j 1 , j 2 ) = (p j1 + p j2 ) 2 . (5.39) For the heavy neutrino masses under consideration, we see that the peaks of the invariant mass spectra occur at M (j 1 , j 2 ) ∼ 1000 GeV − 1200 GeV, with a peak position at larger M (j 1 , j 2 ) for larger m N . A narrow collection of events at M (j 1 , j 2 ) 500 GeV is also observed. These low-mass events are attributed to instances of one forward jet possessing relatively low p T while another jet undergoes a hard q * → qg splitting. In such cases the (qg) pair can be identified as the leading jet pair but still return a small M (j 1 , j 2 ) since this corresponds to the (q * )-system's virtuality, which is favored to be small in massless QCD. The dependence on m N indicates that the hadronic activity is not completely decoupled from the W ± W ± → µ ± µ ± sub-process, and therefore can potentially offer a handle on determining the value of m N . This is relevant given the mild scale uncertainty bands.
The M (j 1 , j 2 ) spectra point to the signal process being driven by valence quark-valence quark scattering involving large momentum fractions, i.e., x B > M (j 1 , j 2 )/ √ s ∼ 0.1. In comparison to the p µ T dis- tributions of figure 3, which show charged lepton momenta reaching a few hundred GeV, we see that comparable momentum fractions are propagated into the W ± W ± → µ ± µ ± sub-process. For example: estimating the incoming W boson energies by those of the muons, E W ∼ E µ ∼ p µ T ∼ 100 GeV − 300 GeV, and the outgoing quark energies from the invariant mass of the two leading jets, which are also back-to-back, E out q ∼ M (j 1 , j 2 )/2 500 GeV − 1000 GeV, then the typical momentum fractions carried by the W reach Moving onto figure 6(d), we show the pseudorapidity difference between the two leading jets, defined as ∆η(j 1 , j 2 ) = η j1 − η j2 . (5.41) We report several notable features. First is the symmetric behavior around ∆η = 0, which stems from having a symmetric beam configuration. Second is that most of the phase space populates the region where |∆η| 2 and appears independent of heavy neutrino masses. Third is the presence of a modest collection of events with |∆η| 1. Such events are consistent with the low-mass distribution in figure 6(c) originating from q * → qg and g * → qq splittings with relatively small opening angles. Beyond one-and two-particle observables are those sensitive to the global activity of the pp → µ ± µ ± jj + X process. In particular, we consider in figure 7(a) the scalar sum of p T over all jets in an event (H T ), and in figure 7(b), the scalar sum of the p T of all recon-structed particle candidates (reco) (X T ), In the first case and for all heavy neutrino masses we observe that H T peaks at H T ∼ 100 GeV and uniformly decreases for larger values of H T . As the net contribution of the two (three) leading jets in the signal scales as 44) the H T distribution suggests the presence of little highp T hadronic activity beyond these leading objects. For the X T case, we observe a slight dependence on N 's mass, with the distributions' maxima occurring at X T ∼ 600 GeV − 750 GeV and tending towards larger values for larger m N . We attribute this sensitivity to the dependence of muon p T , as seen in figure 3, which also peaks at larger values for increasing m N . In comparing X T to the (µ ± µ ± jj(j))-system itself, we find that the scalar sum of H T and the same-sign muon pair p T , Finally, we consider observables that measure the relative amounts of hadronic and leptonic activity in a given event. Such quantities are sensitive to the color structure of hard scattering processes [30,90,132,133], and hence employable in dynamic jet vetoes for color-singlet signal processes. While a full exploration of jet vetoes in W ± W ± → µ ± µ ± is outside our scope, we consider as representative cases in figures 7(c) and 7(d) respectively, the ratio of H T and the leading charged lepton p T : and the ratio of X T and the leading charged lepton p T : In the first case, we observe that a majority of the phase space sits well below r H T µ1 = 1. This indicates that on an event-by-event basis more transverse momentum is carried by the leading muon than in all jets combined. This is unlike diboson and top quark processes where the situation is reversed [30,90,133]. We find that the ratios peak just above r H T µ1 ∼ 0.25 and are largely independent of the heavy neutrino masses that we consider. This is also consistent with naïve estimations from the individual H T and p µ1 T distributions, which suggest In the r X T µ1 case, we observe a sharp cut-off at r X T µ1 ∼ 2 that is largely independent of the heavy neutrino's mass. This can be tied to the disparity of momentum scales between leptons and jets. In particular, since p µ1 T , p µ2 T p j k T , one finds that the ratio scales roughly as That p µ1 T ∼ p µ2 T is consistent with the back-to-back trajectories found in the azimuthal separation in figure 6(a).

VI. SENSITIVITY AT THE LHC AND HL-LHC
In this section we estimate the discovery potential of heavy Majorana neutrinos in same-sign W W scattering at the LHC and its high-luminosity upgrade. After summarizing our simulated detector setup in section VI.1, we build our event selection menu in section VI.2, and present our findings in section VI.3. Our analysis includes signal and background processes that are normalized to an integrated luminosity of L = 300 fb −1 for the √ s = 13 TeV LHC, and to L = 3 ab −1 for the HL-LHC.

VI.1. Detector Modeling and Particle Identification
Particle objects considered throughout the analysis are defined using the ATLAS configuration card available from the Delphes repository. In the results that follow the ATLAS card was modified to construct jet candidates from hadronic activity using the anti-k T sequential clustering algorithm with a distance parameter R = 0.4. This value is more widely used in recent ATLAS data analyses than the default value of R = 0.6. All the other parameters in the configuration card, which include particle identification and mis-tagging efficiencies as well as the fiducial geometry definition, are left unchanged.
As summarized in section III.1, the impact of pileup is assumed to be subtracted from events as one would do with real data. While we neglect the presence of additional low-p T , pileup jets in our samples, the impact on particle resolution is at least partly encapsulated in the momentum smearing routines in Delphes.

VI.2. Event Selection
Our analysis is designed to be as simple as pragmatically possible. We do this to establish a baseline sensitivity and discovery potential at the LHC that broadly covers Majorana masses spanning m N = 50 GeV − 20 TeV. Investigating improvements that target specific mass regimes is left to future work. As discussed in section VII.1, we do not exploit all the kinematic characteristics reported in section V.2. We omit, for example, cuts on E miss T or hadronic activity that are well-established handles in searches for LNV at colliders [6,25,30]. Thus, a more tailored analysis by ATLAS or CMS should yield improvements over the outlook presented here.
To identify our LN-violating collider signature, which is characterized by two same-sign muons and at least two jets, we first apply a loose event selection (called pre-selection in the following). This reduces the number of background processes that must be considered while keeping a high selection efficiency for the signal process.
The pre-selection also includes requirements needed to ensure the near 100% efficiency of the inclusive, singlemuon trigger chains used to record collision data in AT-LAS [134] and CMS [135] during Run 2 and the future Run 3 of the LHC. A summary of pre-selection requirements is listed in the top of Table III. At pre-selection we require that events have exactly two isolated muon candidates. Muon candidates must have the same electric charge Q and reside within the fiducial volume of |η| < 2.7. The leading muon must have p T > 27 GeV in order to ensure full efficiency of the muon trigger, whereas the sub-leading muon p T just needs to be above the reconstruction threshold of p T > 10 GeV. Events with additional lepton candidates, including hadronically decaying τ leptons, are vetoed.
At least two jets with p T > 25 GeV and |η| < 4.5 must be present in each event. The invariant mass of the system constituted by the two highest p T jets passing these criteria, M (j 1 , j 2 ), must also exceed 700 GeV. This suppresses interfering sub-processes (see section IV.1) and enriches the signal sample with a topology that corresponds to the scattering of weak vector bosons.
For the heavy neutrino mass range under consideration, we report that about A ∼ 20% to 40% of signal  (3.3), the preselection-level cross section σ Pre. (third column), and the pre-selection acceptance rate A = σ Pre. /σ Gen. . After pre-selection cuts, we anticipate that the leading background processes consist of mixed EW-QCD production of W ± W ± jj, pure EW production of W ± W ± jj, and the inclusive diboson(+jets) spectrum W ± V + nj with V ∈ {γ * /Z ( * ) }. For compactness we label these: Our MC modeling of these backgrounds is described in section IV.2 and their generator-level rates are summarized in Table I. We have checked that other processes satisfying the same-sign muon signature do not appreciably survive pre-selection. For example: using the Kfactors of Ref. [136], we estimate that the rate for the ttW ± → 2µ ± + X process after pre-selection cuts, but minus the M (j 1 , j 2 ) requirement, is σ ttW ± ∼ 1.2 ab. The M (j 1 , j 2 ) criterion reduces this an order of magnitude. We acknowledge that we do not fully account for backgrounds that are difficult to simulate accurately with MC simulations alone. This includes opposite-sign muon pairs in which one muon is reconstructed with the incorrect charge, or events where one of the same-sign muons has a "fake" origin, e.g., a muon candidate originating from a jet. While such backgrounds are sub-dominant in the dimuon final state, this is less so for other lepton flavors [30,[137][138][139]. Whatever the final state, such backgrounds should be investigated carefully by experiments through dedicated studies based on collision data, as usually done in collider searches for LNV [18,20]. T after pre-selection cuts for the signal process assuming representative heavy neutrino masses mN = 750 GeV (darkest cross), 1.5 TeV (dark cross), and 5 TeV (light cross). Also shown are the inclusive W ± V (light), EW W ± W ± jj (dark), and QCD W ± W ± jj (darkest) backgrounds, as well as their sum (black cross).
To help define our analysis's signal region, we present in figure 8 the distribution of (a) M (j 1 , j 2 ) and (b) p µ2 T for the signal process after pre-selection cuts, assuming representative heavy neutrino masses m N = 750 GeV (darkest cross), 1.5 TeV (dark cross), and 5 TeV (light cross). We also plot after pre-selection cuts the W ± V (light), EW W ± W ± jj (dark), and QCD W ± W ± jj (darkest) backgrounds, as well as their sum (black cross). The curves are normalized to L = 300 fb −1 .
Qualitatively, we observe that background processes tend towards smaller values of transverse momentum and invariant mass while the signal process tends towards larger values and exhibit broader, wider distributions. More quantitatively, we observe in figure 8(a) that background processes peak at M (j 1 , j 2 ) 800 GeV and taper off for larger invariant masses. This contrasts with the signal samples, which peak at M (j 1 , j 2 ) 900 GeV, plateau for a couple hundred GeV, and then gradually fall off. While the lightest heavy neutrinos benchmarks stay above the SM background for most all values of M (j 1 , j 2 ), we observe that the heaviest benchmark mass at m N = 5 TeV only exceeds the background for M (j 1 , j 2 ) 4.5 TeV. Values of active-sterile mixing below unity will naturally worsen this separation power.
In figure 8(b) we observe that all backgrounds peak at p µ2 T ∼ m V /2 ∼ 40 GeV − 45 GeV, and quickly dissipate at higher p T . As anticipated, this shows that backgrounds are driven by resonant weak boson production, though not exclusively. Signal rates become more prominent for p T 50 GeV − 100 GeV. For heavy neutrinos masses beyond a few hundred GeV we report high selection efficiencies when requiring p µ2 T above this range, but less so for lower masses. In this regime, developing an alternative analysis strategy may increase the sensitivity but goes beyond the scope of this work.
For p µ2 T 300 GeV, we find that the total background rate is strongly suppressed. Subsequently, due to its simplicity, we define our signal region by requiring, in addition to pre-selection requirements, that both same-sign leptons carry p T above 300 GeV. We summarize this in the bottom of Table III. For the heavy neutrino masses under consideration, we find that about ε ∼ 15% to 80% of pre-selection signal events survive signal region requirements. For representative masses, we report in last column of table IV the signal rate cross section σ SR and the corresponding selection efficiency ε = σ SR. /σ Pre. . After all selection cuts, we find that the total background rate reaches about σ All cuts b ≈ 2.35 ab. For each background and their sum, we list in Table V the expected number of background events after full selection for the nominal LHC scenario (LHC) with L = 300 fb −1 as well as for the high-luminosity scenario (HL-LHC), where our estimate is computed by simply rescaling the luminosity to L = 3 ab −1 . At the LHC, less than one background event is expected to pass the selection.

VI.3. Results
To quantify the expected excess number of W ± W ± → µ ± µ ± signal events over the number of SM background events, we follow the recommendations of Ref. [140] and employ asymptotic distributions of test statistics. In par- . Expected 95% CL sensitivity at the √ s = 13 TeV LHC (300 fb −1 ) and HL-LHC (3 ab −1 ) on the squared activesterile mixing element |VµN | 2 as a function of heavy neutrino mass following the W ± W ± → µ ± µ ± analysis described in the text. Also shown are the direct limits (35.9 fb −1 ) set by CMS using the CCDY and W γ fusion channels [13], an extrapolation of the CMS to the HL-LHC, and indirect limits [68]. ticular, we define our signal significance Z as [141,142]: log y , with (6.5) , and y = 1 + . (6.6) Here, n = n s + n b is the total number of observed events, n s = L × σ SR s is the number of signal events expected for an integrated luminosity of L and signal region rate σ SR s . The quantity n b = L × σ SR b is the number of background events expected for a signal region rate σ SR b , and δ b is the uncertainty in n b . Based on the background uncertainties recorded in Table I and the backgrounds that we neglect, we conservatively take our uncertainty in the background to be 20%, i.e., we set δ b = 0.2.
As discussed in Ref. [140], the significance estimate Z is consistent with a Poisson-counting likelihood where the background-rate nuisance parameters are constrained by auxiliary Poisson measurements [141]. This constraint is performed, for example, by using control samples enriched with background events. Defining such control samples, which are ultimately employed in the likelihood fits that constrain the normalization of backgrounds in the signal region, is beyond our scope. They are, however, commonly carried out in LHC analyses by choosing control samples in region of phase spaces as close as possible to the signal region, but where no signal is expected.
We estimate the LHC's discovery potential (Z = 5) and exclusion limits at 95% confidence level (CL) (Z = 1.96 in the Gaussian limit, which we approximate with Z = 2) to heavy neutrinos and their active-sterile mixing |V µN | by setting n s to be where n 0 s is the expected number of signal events when |V µN | is fixed to unity. For each mass point, we then set Z = 5 (2) for discovery potential (exclusion sensitivity) and solve the resulting equation for |V µN | 2 . As a function of heavy neutrino mass, we show in Figure 9 and Table VI the discovery and exclusion limits for our luminosity benchmarks at L = 300 fb −1 and 3 ab −1 . In figure 9, we also show the direct limits (35.9 fb −1 ) set by CMS using the CCDY and W γ channels [13], an extrapolation of the CMS to L = 3 ab −1 , and indirect limits [68].
We report that |V µN | 2 0.05 − 0.5 (0.03 − 0.3) can be probed at 95% CL for m N = 50 GeV − 20 TeV with L = 300 fb −1 (3 ab −1 ). Under the assumption that the mass scale of one or more heavy neutrinos is much heavier than collision scales at √ s = 13 TeV, then in analogy to interpretations of searches for 0νββ decay, the LHC expected sensitivity at 95% CL can be expressed as Due to higher background rates, we anticipate slightly worse sensitivity for the e ± µ ± jj and e ± e ± jj channels. For final-states involving hadronic taus, we anticipate even weaker (but still comparable) sensitivity due to tagging efficiencies. Dedicated studies of these complementary signatures are strongly encouraged.
We find that the proposed analysis has a strong potential to significantly extend the current sensitivity of ATLAS and CMS searches for resonant heavy neutrino masses beyond m N ∼ 750 GeV and up to masses at the O(10 TeV) scale, as shown in figure 9. Similar to the existing ATLAS and CMS searches, the proposed analysis however does not reach a sensitivity comparable to the limits set by indirect precision measurements (see section II.2). The analysis does, however, offer a direct test of the 0νββ decay mechanism in lepton flavor configurations that are not accessible at nuclear energy scales.

VII. OUTLOOK
Discovering the W ± W ± → ± i ± j process in LHC collisions would present unambiguous evidence that LN is violated at the TeV scale, and have far-reaching repercussions for both theory and experiment. In light of the encouraging sensitivity reported in section VI.3, we now briefly consider the outlook for further improvements to our proposed experimental analysis (section VII.1), as well as the possible application of our work to other scenarios, including when LN is conserved (section VII.2).

VII.1. Improving the Experimental Analysis
The analysis cuts chosen and outlined in section VI.2 include only a simple set of selection criteria, which were derived to obtain a good significance for a large range of heavy neutrino masses and with the aim of being robust against the effects of finite detector resolution. Obviously, the selection can be optimized for individual neutrino masses. This is especially true for the lower mass range, where our proposed analysis can add sensitivity to current searches for resonant heavy neutrinos. Such improvements can be roughly grouped into those which further suppress the W ± W ± or diboson backgrounds.
An optimized analysis that also takes into account the resolution and acceptance of specific sub-detectors should considerably improve our sensitivity estimate. An obvious choice would be to use E miss T in the event selection. E miss T would especially help reduce the impact of W ± W ± production, which has otherwise the same topology as the signal, but was omitted since it is particularly sensitive to detector resolution. As examined in section V.2, the heavy neutrino signal is characterized by high muon momenta that can reach the TeV scale. Even for a small fractional mis-measurement of muon momentum this can induce a considerable amount of E miss T . Hence correlating the E miss T and muon momentum (or similarly jet momen-tum) in the selection would be a way to mitigate some of the resolution effects.
W W scattering processes, such as the heavy neutrino signal considered in this paper, commonly feature suppressed central hadronic activity. As a consequence, implementing static jet vetoes [108][109][110][111], or their dynamic counter part [30,90,133,143,144], can be exploited to further reduce diboson and top quark processes (or in general all non-VBF backgrounds). One can optimize the corresponding selections based on detector efficiency and resolution for low momentum jets as well as take into account LHC pileup conditions. Since the signal process discussed in this paper is modeled at NLO in QCD with parton shower-matching, it can be used to study improvements connected to central hadronic activity.
Lastly, our projections for the HL-LHC consisted of only a re-scaling of the results obtained for a smaller data set. However, the planned detector upgrades for ATLAS [145] and CMS [146] will also allow one to refine the analysis's selection criteria. In particular, the extended coverage of the tracking system will be highly relevant for W W scattering processes due to the use of tracking information for jets in the forward region.

VII.2. Applications to Other Seesaw Searches
In this study we have focused exclusively on the LNviolating process pp → ± i ± j jj, when mediated by samesign W W scattering and the t-channel exchange of a heavy Majorana neutrino, as shown in figure 1. In a complementary fashion, the LN-conserving process pp → ± i ∓ j jj, which proceeds through opposite-sign W W scattering, is also possible [26,36]. One could anticipate that the differences in helicity inversion (see section V.1) between the LN-violating W ± W ± → ± i ± j sub-process and the LN-conserving W + W − → ± i ∓ j sub-process results in substantially different kinematic distributions. However, for heavy neutrinos in the range of m N = 750 GeV to 5 TeV, this may not be the case.
As reported in section V.2, we found that the W ± W ± → ± i ± j sub-process in pp → ± i ± j jj behaves like a factorizable system with kinematics that are nearly independent of the hadronic environment. This means that much of the charged lepton kinematics are driven more by momentum conservation in 2 → 2 scattering than some complex spin-correlation. It is arguable that many of the kinematic leverages over backgrounds that we find, e.g., large dijet invariant masses and back-toback charged lepton trajectories, will also hold for the LN-conserving case. While the pp → ± i ∓ j jj collider signature inherently has a much larger background rate than the LN-violating one, we stress that lepton flavor violation is forbidden in the SM. Therefore, requiring that i = j and that E miss T is small, as done for example in Ref. [147] for low-scale Type I Seesaw models [50][51][52][53][54][55], can significantly reduce SM backgrounds.
As a final remark, we comment on the applicability of our analysis to other neutrino mass models. This include, for example, the Type II Seesaw model [49,[148][149][150][151], wherein the LN-violating, VBF sub-process W ± W ± → ∆ ±±( * ) → ± i ± j can occur through a possibly non-resonant, s-channel, doubly charged Higgs boson ∆ ±± , as well as the Left-Right (LR) Symmetry model [152][153][154][155][156], wherein the LN-violating, VBF subprocess W ± R W ± R → ± i ± j can proceed through righthanded W ± R gauge bosons and a t-channel Majorana neutrino. For the LR scenario this is interesting in the event that resonant production of W R is not within the kinematic reach of the LHC [157,158]. As both scenarios can mimic our pp → ± i ± j jj collider signature, its discovery does not automatically prove the existence of RH neutrinos. On the other hand, as both processes occur through the scattering of two color-singlet, massive gauge bosons, most of the color and Lorentz structure remains identical to the original case that we study. Therefore, we anticipate again that the VBF sub-processes approximately factorize, resulting in charged lepton kinematics that resemble those presented in section V.2. As a result, the collider analysis that we propose can readily and justifiably be recast in terms of the Type II and LR models.

VIII. CONCLUSIONS
Motivated by the possible non-conservation of LN in nature, we have investigated the potential to search for heavy Majorana neutrinos in same-sign W ± W ± scattering at the LHC and the HL-LHC. The experimental signature of two forward jets from VBF, two same-sign leptons, and the lack of substantial missing transverse momenta is interesting in its own right as, to our knowledge, it was not yet explored experimentally at the LHC.
As a benchmark scenario we use the Phenomenological Type I Seesaw model with two key aspects to be probed experimentally: the mass m N of a mostly sterile neutrino N and its mixing with the active neutrinos |V N |. Current searches at the LHC target resonant production modes, such as the Drell-Yan and W γ fusion mechanisms, which have the advantage of a factor |V N | 2 less suppression compared to the W ± W ± → ± i ± j channel. They however suffer from rapidly falling scattering rates at increasing heavy neutrino masses due to matrix element and phase space suppression. For these reasons LHC searches that employ resonant production modes can only probe masses up to m N = 3 − 4 TeV [30].
To conduct this study, we developed in section IV simulation prescriptions at NLO in QCD with parton showermatching for both the VBF signal process and backgrounds based on the HeavyN UFO libraries and the MG5aMC simulation suite. We then extensively studied in section V the phenomenology of the W ± W ± → ± i ± j process at the amplitude and differential levels. We find that "bare" cross section for the full, 2 → 4 signal process at NLO peaks for heavy neutrino masses of around 1 TeV and can reach up to σ/|V N | 4 ∼ 10 fb at √ s = 13 TeV.
Apart from the large rapidity gap between the two leading jets and large dijet invariant mass, the signal also features very high lepton momenta, among other characteristics, which can be exploited for an effective background suppression.
In section VI we designed our collider analysis, employing the Delphes framework to simulate the response of a typical LHC detector. Our analysis was deliberately kept simple and considers only final states with same-sign muon pairs to obtain reliable and robust results. Accordingly, dedicated analyses exploiting the suppressed QCD radiation in VBF processes, the angular separation of the same-sign lepton pair, or the correlation between the measured missing transverse momenta and very high-p T leptons should improve on our projected sensitivity.
In section VI.3 we show that with the LHC Run 2 and expected Run 3 data sets, |V N | 2 can be probed down to 0.1-0.5 at 95% CL for heavy neutrino masses in the range m N = 1 − 10 TeV, and that masses at m N = 20 TeV can be probed for |V N | 2 down to 0.8. At the HL-LHC, this can be improved by a factor of 2 to 5. We find that the W ± W ± fusion channel extends significantly the current mass reach based on resonant production modes and adds valuable sensitivity to the masses above a few hundred GeV. Finally in section VII, we give an outlook on areas where our proposed analysis can be improved as well as on complementary applications of our results. In this appendix we derive constraints on heavy Majorana neutrinos that arise from direct searches for nuclear 0νββ decay as reported in equation (2.7). To do this, we assume that the decay is solely mediated by the exchange of heavy states N k of mass m N k that couple to SM particles according to the Lagrangian of equation (2.3). We work in the standard factorization picture [159][160][161][162]. This stipulates that the transition rate for the decay process of nucleus (A, Z) into nucleus (A, Z + 2), can be expressed as a product of the two-body phase space factor G 0ν for the (e − e − )-system; a nuclear matrix element (NME) A; and the propagators for the states N k . Under the assumption that the 0νββ decay process is mediated only by the t-channel exchange of W bosons and N k , the NME simplifies to A ≈ A N . Explicitly, these are related to the decay half-life (T 0ν 1/2 ) by The proton mass m p ≈ 0.938 GeV is introduced to render A N dimensionless. Typical momentum transfers in nuclear 0νββ decay are of the order |t| ∼ O(0.1 GeV).
As we are interested in EW-scale and TeV-scale neutrinos, this allows us to expand equation (A2) and obtain By neglecting O |t|/m 2 N k contributions we can invert equation (A2) and translate an experimental lower bound on T 0ν 1/2 into an upper bound on the mixing-over-mass ratio of Majorana neutrinos. In the following we focus on the 76 Ge→ 76 Se+2e − transition, and consider the recent experimental limits by the GERDA experiment [67,163].
Following Ref. [164], we use the NMEs of Refs. [165,166]. These employ the so-called Self-consistent Renormalized Quasiparticle Random Phase Approximation (SRQRPA), and make use of two potential models to describe the nucleon-nucleon interactions, namely the Argonne and Charge Dependent Bonn (CD-Bonn) models. The calculations moreover rely on intermediate and large single-particle spacing, i.e., eigenstate multiplicity, an axial-vector coupling constant g A = 1.25, and a nuclear radius R = 1.1 fm × A 1/3 . For 76 Ge, we list in Table VII  the A N for the four nuclear potential configurations. The uncertainty in A N is estimated by considering the envelope spanned by the configurations. We use the phase space factor G 0ν = G (0) 0ν g 4 A , as derived by Ref. [167], assuming an axial-vector cutoff of g A = 1.25. While the polarization-dependent component G (1) 0ν is non-zero, its impact on the total 0νββ decay rate vanishes after phase space integration. The derivation of the energy-dependent component G 0 0ν ∝ 1/R 2 in Ref. [167] uses a nuclear radius of R = 1.2 fm × A 1/3 . Hence, for consistency with our NMEs, we rescale it by For 76 Ge we obtain the phase space factor, G 0ν ≈ 6.866 × 10 −15 yrs −1 . (A6) When added in quadrature, the total estimated uncertainty in this number spans about δG 0ν ≈ 7%−9% [167]. This is considerably smaller than the NME uncertainty, and therefore is neglected in our final constraints on N k . After an exposure of E = 127.2 kg-yrs, GERDA reports a lower limit on the 0νββ decay half-life in 76 Ge of [67] T 0ν 1/2 > T GERDA 90% CL = 1.8 × 10 26 yrs at 90% CL.
For the range of NMEs, this translates into the following upper limit on Majorana neutrino masses and mixing: For each NME that we consider we report in Table VII the corresponding limit on the mixing-to-mass ratio. For related discussions on 0νββ decay, see Refs. [168][169][170].
As described in section IV.2.3, our baseline modeling of the diboson spectrum pp → 3 ν + X, provides limited MC statistics when the two same-sign leptons carry p T 100 − 150 GeV but the odd-sign lepton is much softer. To populate this phase space region, we introduce into MG5aMC's phase space integration routines tailored generator-level cuts (p −cut T ) on the p T of the two leading charged leptons, independent of charge. This is in addition to the baseline cuts of equations (4.3) and (4.4); no further cut is applied to the trailing charged lepton.
High-statistics samples with p −cut T = 75 GeV and 200 GeV are stitched to the baseline FxFx1j sample through cuts on the truth-level p T of the sub-leading charged lepton. Within statistical uncertainty, we report that the shape and normalization of the high-p T tails for the leading charged leptons in the stitched sample reproduce those in the baseline FxFx1j sample. In practice, the magnitude of the charged leptons' transverse momenta are first added, then the largest and smallest p T are subtracted to extract the p T of the subleading charged lepton. If either the leading or subleading p T are below the p −cut T threshold (pTlXCut), then the phase space point is rejected.
To ameliorate inefficient phase space sampling associated with our cuts, we increment the boundary of the PDF convolution integral τ min =ŝ/s, where This is inserted just after the enddo closure tag at about L421 and just before the line stot = 4d0*ebeam(1)*ebeam (2) For phase space cuts beyond p −cut T ∼ 150 GeV, we observe a severe instability in phase space integration. As documented elsewhere (see footnote 4), this failure is attributed to inefficient phase space sampling for nonresonant diagrams with massive τ leptons. Hence, for the p −cut T = 200 GeV sample, we import into MG5aMC the model file loop sm-no tau mass, which assumes a massless τ lepton. For looser p −cut T , we find that this results in sub-percent differences in the cross section normalization from the loop sm model file.
formalism, which is akin to collinear factorization in perturbative QCD, the W boson is treated as a parton of the proton. This allows us to focus on the 2 → 2 subprocess where p and λ denote the 4-momenta and helicities of external particles. The amplitudes reported here supplement the analytic results reported in section V. They are also complementary to the numerical results reported throughout the main text, which evaluate precisely the full 2 → 4 helicity amplitudes for the pp → ± ± jj + X process using the HeavyN NLO UFO libraries [29] in conjunction with the MG5aMC [75,76] MC event generator (see sections III and IV for related details).
For the above process, we work in the hard-scattering frame, which is equivalent to the W W scattering frame. In this frame, we align coordinate axes such that Here M 2 W W = (p W 1 + p W 2 ) 2 is the invariant mass of the (W W )-system and the remaining invariants are given by Working in the unitary gauge and assuming a clockwise fermion flow of leptons [113,114], the helicity amplitudes in the HELAS basis [115] are given by where the (t ↔ u) term accounts for final-state lepton exchange, and the LN-violating tensor current T µν is u(p 1 , λ 1 )γ µ γ ν P L v(p 2 , λ 2 ) .
We assume the exchange of a single sterile neutrino mass eigenstate N with momentum p N = (p W 1 − p 1 ) and mass m N . For the more general case of multiple heavy neutrinos N k with masses m N k , one would substitute: to capture the interference of multiple propagating mass eigenstates. Similarly, for two final-state lepton flavors, one substitutes, V N k V N k → V 1Nk V 2Nk . Importantly, the spinor and Lorentz index contractions are not modified and therefore are the same for any number of tchannel Majorana neutrino exchanges. Explicit evaluation (and inspection) of T µν shows that the tensor is non-vanishing only when both final-state antileptons carry right-handed polarizations, (λ 1 , λ 2 ) = (R, R). Beyond this, the full matrix element also vanishes identically when both incoming W bosons carry opposite transverse polarizations, i.e., when (λ W 1 , λ W 2 ) = (±1, ∓1), which follows from an orthogonality of d-and p-wave states. For both t-channel (second column) and u-channel (third column) exchanges of a heavy neutrino, we list in table VIII the exact helicity amplitudes as a function of external particle helicity (first column).
We find that t-and u-channel tensor structures for the (λ W 1 , λ W 2 ) = (0, ±1) and (±1, 0) configurations differ simply by a global factor of −1. The t-and u-channel structures for (λ W 1 , λ W 2 ) = (±, ±1) differ by exchanges of sine and cosine functions, whereas for (λ W 1 , λ W 2 ) = (0, 0) there is a minor difference in the polar angle dependence. For a fixed set of helicity polarizations, the matrix element of equation (C8) is obtained by the standard coherent summation of t-and u-channel terms.

Low-mass limit
When the heavy neutrino and W boson masses are both small compared to the (W W )-scattering scale, i.e., m N , m W M W W , we can expand each of the squared matrix elements in powers of the ratios Doing so reveals that remaining transverse-transverse permutations, i.e., (λ W 1 , λ W 2 ) = (±1, ±1), as well as LHlongitudinal channels, (λ W 1 , λ W 2 ) = (−1, 0) and (0, −1), either vanish or are sub-leading. In the latter cases, we see the emergence of a helicity suppression that can compete or overcome longitudinal polarization enhancements, which scale as ε µ (k, 0) ∼ k µ /m W + O(m W /k 0 ).
the alignment of the (W W )-system's angular momentum with that of the dilepton system and a single longitudinal polarization enhancement. Explicitly, we obtain This expression does not account for the (1/2!) multiplicity factor for identical final-state particles. We find that the leading polarization configuration in this kinematic limit is the longitudinal-longitudinal channel, (λ W 1 , λ W 2 ) = (0, 0). We attribute its survival in the expansion to the double longitudinal enhancement and is To build the partonic, 2 → 2 cross sectionσ we employ the standard relationship between scattering rates and matrix elements. This is given by the phase space integral where the 2-body phase space volume measure is dP S 2 (p W 1 + p W 2 ; p 1 , p 2 ) = d cos θ 1 dφ 1 2(4π) 2 β W , (C17) and the differential scattering rate is The symmetry factor S = 3 2 · 2! accounts for spinaveraging over initial-state W polarizations and identical, final-state particles, while the incoherent summation runs over all external helicities {λ W 1 , λ W 2 , λ 1 , λ 2 }. The velocity factor β W = √ 1 − 4r W accounts for the masses of incoming beam particles.
After phase integration over the azimuthal direction, the leading contribution to the polar distribution of + 1 in W + W + → + + scattering is given analytically by After integration over the polar angle, the total rate iŝ For this limit, we report agreement between this expression and numerical evaluations of the same 2 → 2 process using the HeavyN model with MG5aMC. For multiple heavy neutrinos coupling to potentially different charged lepton flavors, the above generalizes tô

High-mass limit
When the heavy neutrino mass is large compared to the W ± W ± → + + scattering scale, i.e., when m N M W W , m W , one can work in the decoupling limit [120] and treat the exchange of N as a point-like, contact interaction. Formally, this entails expanding the heavy neutrino propagator such that Inserting this expansion into the amplitudes listed in table VIII reveals a strong destructive interference among most of the helicity permutations. In particular, the only non-vanishing channels correspond to those with incoming W bosons carrying identical polarizations, i.e. (λ W 1 , λ W 2 ) = (0, 0) and (±1, ±1). Explicitly, the matrix elements for these configurations are given by In comparing the two expressions one can see the impact of the longitudinal polarization enhancements, which are responsible for the relative factor of (M W W /m W ) 2 . After squaring and integrating over the azimuthal angle, we obtain the leading contributions in the decoupling limit to the polar distribution of + 1 in W + W + → + + scattering. For each polarization channel, this is: We report good agreement between these expressions and numerical evaluations of helicity-polarized cross sections in this kinematic limit using the HeavyN model with MG5aMC in conjunction with the formalism of Ref. [171].
After integrating over the polar angle and averaging over all W boson helicities, the total scattering rate iŝ For n R heavy neutrinos coupling to potentially different charged lepton flavors, the above generalizes tô