Heavy neutral leptons in effective field theory and the high-luminosity LHC

Heavy neutral leptons (HNLs) with masses around the electroweak scale are expected to be rather long-lived particles, as a result of the observed smallness of the active neutrino masses. In this work, we study long-lived HNLs in $N_R$SMEFT, a Standard Model (SM) extension with singlet fermions to which we add non-renormalizable operators up to dimension-6. Operators which contain two HNLs can lead to a sizable enhancement of the production cross sections, compared to the minimal case where HNLs are produced only via their mixing with the SM neutrinos. We calculate the expected sensitivities for the ATLAS detector and the future far-detector experiments: AL3X, ANUBIS, CODEX-b, FASER, MATHUSLA, and MoEDAL-MAPP in this setup. The sensitive ranges of the HNL mass and of the active-heavy mixing angle are much larger than those in the minimal case. We study both, Dirac and Majorana, HNLs and discuss how the two cases actually differ phenomenologically, for HNL masses above roughly 100 GeV.

So far, searches for physics beyond the standard model (BSM) at the Large Hadron Collider (LHC) have focused on promptly decaying heavy new states. However, no concrete discovery of any such field has been announced yet, despite the roughly 140 fb −1 of statistics accumulated at both CMS and ATLAS. Perhaps mostly for this reason, attention in recent years -from both experimental and theoretical fronts -has shifted towards the study of fields that are light and very weakly coupled to SM fermions, such that they have a relatively long lifetime [1][2][3].
Long-lived particles (LLPs), once produced at a collider, can lead to variety of exotic signatures as their decay products are displaced from their production position [1]. For instance, charged particles nearly mass degenerate with a neutral, stable dark matter candidate can lead to a disappearing track signature. An example of such exotic signature is predicted in the Minimal Supersymmetric Standard Model (MSSM) with a compressed electro-weakino spectrum (for searches at the LHC, see Refs. [4][5][6]). Electrically neutral LLPs, on the other hand, can lead to signals with displaced vertices. The literature abounds with candidate models. For example, in R-parity-violating supersymmetry [7][8][9] the lightest SUSY particle can be long-lived [10,11] because of the smallness of the R-parity-violating coupling due to the small neutrino masses.
A particularly simple and popular class of neutral LLP models are the so-called 'portal' models. In this class of models, a new scalar, pseudoscalar, fermion or vector field is introduced, connecting the SM sector with a hidden (or dark) sector. In particular, the fermionic-portal mediator is conventionally called a 'heavy neutral lepton' (HNL) in the context of LLPs. Event numbers in such models are expected to be very small. Nevertheless, LLP searches benefit from considerably lower backgrounds than prompt searches. In addition, given that the LHC will deliver up to 3 ab −1 of statistics over the next decade and a half, large parts of untested parameter space of LLP models could be probed.
In the minimal HNL scenario, there is a small mixing of unspecified origin of the HNLs with the active neutrinos. This setup is motivated by the neutrino masses (and mixings) observed in oscillation experiments [12] (see also Refs. [13,14]), and the HNLs could be identified with the sterile neutrinos in a 'seesaw mechanism' [15][16][17][18][19]. In the classical seesaw type-I, for O(1) Yukawa couplings sterile neutrinos are expected to be as massive as The best systematic way to study such "non-minimal" HNLs is to apply effective field theory (EFT). In EFT, non-renormalizable operators (NROs) with mass dimension d > 4 respecting the gauge symmetries of the SM are introduced on top of the standard renormalizable terms. NROs are suppressed by a new physics (NP) scale Λ that is assumed to be higher than the electroweak scale. In this work, we denote the EFT of the SM augmented with sterile right-handed (RH) neutrinos, N R , as N R SMEFT [60][61][62] (see also Refs. [63][64][65]). A complete classification of d ≤ 7 operators containing N R can be found in Ref. [62]. 1 At d = 5, there are two operators involving N R . Both of them violate lepton number. Their phenomenology, and in particular, that of RH neutrino magnetic moments, was originally investigated in Ref. [61]. The effects of the d = 5 Higgs-N R operator leading to a decay of the Higgs boson to two HNLs (each of which subsequently decays to a charged lepton and two jets via active-heavy neutrino mixing) were studied in Ref. [67].
Cross sections for lepton number conserving operators with d = 6 were initially discussed in Ref. [60] in the context of Tevatron and LHC. These operators have been further analysed in Ref. [68] in the context of both low-energy and collider experiments. In the last few years the phenomenology of N R SMEFT was studied in different regimes. In particular, in Ref. [69], constraints on the four-fermion operators involving N R have been set from LHC searches for associated production of a charged lepton with missing energy, monojet searches, as well as pion and tau decays, assuming N R to be stable at collider scales. A different regime, in which N R 's decay promptly to an active neutrino and a photon via a dipole operator 2 , was studied in Ref. [72] for operators involving the Higgs boson, and in Ref. [73] for the four-fermion operators. Further collider studies of N R SMEFT include Refs. [74][75][76][77][78].
On the low-energy front, Ref. [79] has derived constraints on the d = 6 operators coming from charged lepton flavor violation (CLFV), beta decays and coherent elastic neutrinonucleus scattering (CEνNS). The constraints originating from the CEνNS as well as from meson decays have been also obtained in Refs. [76,80,81] and those from CLFV in Ref. [82]. Furthermore, the effects of non-renormalizable interactions on neutrinoless double beta decay were investigated in Ref. [83]. Very recently, Ref. [84] has shown that, under certain conditions, a d = 6 operator inducing RH leptonic charged current can account for the observed values of the muon and electron anomalous magnetic moments. Finally, the one-loop renormalization of d = 6 operators containing N R was performed in Refs. [85][86][87] (see also Ref. [63] for earlier partial results).
Most closely related to this work, Ref. [53] has studied long-lived HNLs lighter than about 5 GeV, in the framework of N R SMEFT, taking advantage of the copious production of bottom and charm mesons at the LHC. More concretely, in each of the various theoretical scenarios studied in that work, two separate EFT operators each with one single HNL are switched on together with the underlying electroweak charged and neutral currents, 1 For a basis of the d = 6 operators containing N R and independent off shell, i.e. before applying equations of motion, called also Green basis, see Ref. [66]. 2 For a study of different decay modes of a Majorana HNL triggered by effective interactions see Refs. [70,71]. mediating the production and decay of the HNLs. Sensitivity reaches have been worked out for the fixed-target experiment SHiP [88,89], the local ATLAS detector [90], as well as the far detectors.
In the present work, we focus on operators containing a pair of HNLs and a pair of quarks. Such interactions can enhance the production of the HNLs from parton collisions with respect to the minimal HNL scenario. In these scenarios, the decay of HNLs is mediated largely dominantly by the usual mixing angle 3 , allowing to probe HNL masses well beyond the B-meson threshold. In this theoretical framework, we derive the sensitivity reach of both ATLAS and possible far detectors for the NP scale Λ, the active-heavy mixing angle, as well as the HNL mass. As we will show, Λ in excess of 10 TeV and mixing angles below even type-I seesaw expectations can be probed in this setup.
The rest of this paper is organized as follows. In section II, we introduce the model basics for HNLs including both N R SMEFT and discussing briefly some ultraviolet (UV) complete models. In section III, we describe in detail the simulation procedure and the analysis we perform for both ATLAS and far-detector experiments. Our numerical results are presented in section IV. Finally, in section V, we summarize the work and provide an outlook. We added two appendices. In the first, we present the effects of a possible timing cut at ANUBIS on our analysis, while the second discusses the reinterpretation of a CMS search for HNLs [27] within our theoretical setup.

II. EFFECTIVE FIELD THEORY AND HEAVY NEUTRAL LEPTONS
In this section, we present some basic theoretical aspects of our setup. First, as a motivation, we briefly discuss different variants of electroweak seesaw mechanisms. We then turn to d = 6 NROs in "N R SMEFT", i.e. the EFT of the SM extended by (three generations of) HNLs. Next, we briefly discuss some UV completions which can generate the d = 6 operators in the effective theory. While we will present our numerical results only in terms of EFT parameters, understanding the UV models has been a necessity for our calculation, because of some limitation in MadGraph5 [91][92][93], as we will explain in section III.
A. Heavy neutral leptons, neutrino mixing, and neutrino mass Adding HNLs to the SM is usually motivated as an explanation for the observed (mostly) active neutrino masses. The simplest (and most well-known) model realization of this idea is the type-I seesaw mechanism [15][16][17][18][19]. In this case, the sterile neutrinos are Majorana particles. However, there are other realizations of the seesaw mechanism, such as the inverse [20] or linear seesaw [21,22]. In these models, the heavy sterile neutrinos form quasi-Dirac pairs [94][95][96]. Since there are some interesting differences in our numerical results for Dirac and Majorana cases (see section IV), here we want to briefly discuss how either of the two possibilities can be realized in concrete models.
If we add to the SM particle content two sets of singlets, call them N R and S L , the most general mass matrix of the neutral fermions, in the basis (ν L , N c R , S L ) 4 , can be written as: This setup contains a number of interesting limits. First of all, if (M R , M M ) µ the fields S L decouple, and we find an ordinary seesaw type-I with three heavy and three light states. Light and heavy neutrino masses as well as active-heavy mixing, in leading order in m D /M M , are given by: where In this simplest possible setup we can assume M M is diagonal,M M , without any loss of generality. All six neutral states are Majorana particles in seesaw type-I. m ν is diagonalized by a (3,3) matrix U , which corresponds to the mixing matrix observed in neutrino oscillations. By V we denote the (3,3) matrix connecting the active and heavy states. Using Eqs. (2) and (4), the mixing between the active and heavy sectors can be trivially estimated to be order of V 2 ∼ m ν /M M . Note, however, that this estimate assumes that there are no cancellations among the different terms contributing to m ν , and if this assumption is dropped, larger values of V 2 can be obtained.
In the limit M M = = 0 and µ M R we find the inverse seesaw, while the case M M = µ = 0 and M R corresponds to the linear seesaw. Consider, for example, the case of the inverse seesaw. Again, in leading order in m D /M R , light and heavy sectors mass 4 Here N c R ≡ CN R T represents the charge conjugate of N R , with C being the charge conjugation matrix.
The number of singlets is, of course, arbitrary. For concreteness, here we consider the "symmetric" situation and add one flavour of N R and S L per SM generation of fermions. All submatrices are then (3,3) matrices. In the numerical part of the work we will consider only one generation, for simplicity. 5 This setup [97,98] is sometimes called "double seesaw" [99] in the literature. matrices and active-heavy mixing are given as: Here, {a, b} is the commutator of a, b, and the calculation is done in the basis where M R is diagonal,M R . In this basis both, m D and µ, can have off-diagonal elements in general. The essence of the inverse seesaw is that, since µ is small compared to M R , V is actually expected to be larger than in the case of the standard type-I seesaw, as now we have V 2 ∼ m ν /µ. More important for us is that in the limit µ → 0, the three pairs of N R and S L form each a Dirac particle. Of course, strictly speaking, µ = 0 is not possible for an inverse seesaw model with a consistent phenomenology, since light neutrino masses vanish in this limit. However, as was shown in Ref. [100], as long as the mass splitting of the heavy states (proportional to µ) is smaller than their decay width, the heavy states are quasi-Dirac neutrinos. This implies that (i) they have only lepton-number-conserving decays, and (ii) their production cross section equals the one for Dirac neutrinos. Note that, while the formulas for neutrino mass and for the mass splitting of the heavy states are different, also in the linear seesaw the heavy neutrinos form quasi-Dirac neutrinos, so the phenomenology is qualitatively similar as for the inverse seesaw. In our model implementations (to be detailed below), the calculation of production cross sections with MadGraph5 considers both possibilities: HNLs as Dirac or Majorana states. Parameter fits to neutrino data [12] (see also Refs. [13,14]) could be easily done for either case: for the type-I seesaw with the well-known Casas-Ibarra parametrization [101] and for more general cases with the results of Refs. [102,103]. However, since in our numerical scans we consider only one HNL within the experimentally accessible window, we will not fit neutrino oscillation data explicitly. Rather, we treat the elements of the active-heavy mixing matrix V as free parameters.
Concerning decays of HNLs to SM particles via mixing, these have been calculated several times in the literature. We will therefore not discuss the details here. In our numerical work we use the decay formulas of Refs. [104,105].
We use S L here as the symbol for the left-handed (LH) partner of N R , motivated by the inverse seesaw discussed above, i.e. the pair (S L , N R ) forms a Dirac particle. We want to  distinguish the LH and RH components of this particle, because in the following we will write down only NROs involving N R . This way, the Dirac case of the NROs resembles more closely the Majorana case. 6 At d = 4 one finds the standard Yukawa terms of the different seesaw realizations. NROs in N R SMEFT were studied previously in the literature at d = 5 [60,61,64], d = 6 [60,63,65], d = 7 [62,106], and very recently at d ≤ 9 [107]. A complete list of operators up to d = 7 can be found, e.g. in Ref. [62]. As already mentioned in section I, at the d = 5 level for Majorana N R there are two new types of operators, in addition to the Weinberg operator [108].
At d = 6 level, one finds four-fermion operators and operators containing two fermions plus either derivatives, Higgses or field strength tensors, see Refs. [60,62]. The latter are not interesting to us here, since they are of minor importance for the production cross sections of N R at the LHC. We therefore do not show these operators in table I, where we list baryon and lepton number conserving d = 6 four-fermion operators, containing either one or two N R . We give the particle contents of the operators, without specifying the family indices, which are in general four index tensors. For each operator structure, we provide the number of independent real parameters for n N = 1 and n N = 3 generations of N R (and n f = 3 generations of the SM fermions). We have checked these numbers using the Sym2Int package [109,110]. Finally, stands for the totally antisymmetric tensor connecting SU (2) L indices.
We note that for Majorana neutrinos there is also a lepton-number-violating operator 6 A complete EFT of the SM extended with a Dirac HNL would contain operators involving N R and/or S L .
We would like to emphasise that in this work we assume the HNL itself to be either Dirac or Majorana in nature, whereas usually in the literature on the N R SMEFT, the Dirac case corresponds to the situation where N R is the RH counterpart of the SM LH neutrino ν L , such that the pair (ν L , N R ) forms a Dirac particle. 7 Further, there are two types of operators which violate both, baryon number and lepton number. We do not use them, and thus, do not list them here.
Not all of the operators listed in table I are equally important for LHC phenomenology. From the six types of operators involving pairs of N R , only the first three will significantly affect the production cross sections for N R at the LHC. 8 It is important to note that operators with pairs of N R can not lead to decays of the lightest N R . This has important consequences for phenomenology, since large production cross sections do not imply fast decay for the lightest HNL, if one considers only these operators.
Operators in the right part of table I contain only one N R . They can lead to sizeable production cross sections for N R too, but also mediate decays of (also the lightest) N R to SM fermions at the same time. Thus, for these operators, cross sections and decay lengths are related, which leads to important upper limits on the cross sections, if we wish to study HNLs as long-lived particles.
Thus, in this work we focus on the N R pair operators with quarks. The corresponding d = 6 Lagrangian reads 9 where a sum over the flavour indices i, j = 1, 2, 3 of quarks and k, l = 1, . . . , n N of HNLs is implied. In the numerical part of the work, for simplicity we consider only one generation of the HNLs and only the first generation of the quarks, and thus, simplify the notation to c 11 dN , c 11 uN , c 11 QN and O 11 dN , O 11 uN , O 11 QN , with "11" referring to the first quark generation.

C. Ultraviolet completions for d = 6 operators
The operators of the effective theory in table I can be generated in the UV by either a Z or via leptoquarks (LQs). It is important to notice that no single BSM particle can generate all operators. Since we are mostly interested in studying the EFT limit, here we will only discuss the Z option and three particular scalar LQs. For a complete list of all possible LQ states for generating d = 6 four-fermion operators see Ref. [79].
Let us start with the LQ case, see table II. These three LQs have been chosen, since they can generate the three different N R pair operators with quarks in the left part of table I. Table II gives the naming convention for the LQs and their quantum numbers under the SM gauge group in the order SU (3) C × SU (2) L × U (1) Y , and indicates which N R pair operator is generated by each of these states in the EFT. The different LQs can have the following Yukawa interactions: Here, we have not written down couplings of the LQs to quark pairs, since these would lead to baryon-number-violating processes, if they were present at the same time as the terms in Eqs. (11)- (13). We note that all the couplings are matrices in flavour space. The first terms in each line of these equations will generate the N R pair operators, while the simultaneous presence of two (or more) terms in each line will generate single N R operators in the EFT. For example, , where again we have suppressed flavour indices; the factor of 2 arises from a Fierz identity. It also generates O LdQN with c LdQN /Λ 2 = g QN g dL /m 2 S Q . Since these two operators depend on different combinations of couplings in the LQ model, it is clear that we can treat them as independent operators. Similar comments apply to S d . Interestingly, S u can generate only a N R pair operator O uN .
The N R pair operators of interest can also be generated by integrating out a heavy Z boson. Here, we do not specify a concrete type of Z , but rather treat it as a massive vector boson with the following interaction Lagrangian: Switching on pairs of couplings (g N , Let us note that in this scenario, the pair production of N R would proceed via s-channel diagrams with a resonance enhancement for m Z = 2m N . Thus, for the EFT description to be valid, m Z should be significantly larger than 2m N , contrary to the LQ case, where the N R production occurs via t-channel diagrams.

A. Model implementations
In order to estimate the sensitivity reach of the various experiments considered in this work, we make use of Monte-Carlo (MC) techniques to perform numerical simulations. The first step, common to all experiments, is the implementation of the different model setups in MadGraph5 [91][92][93]. We built our models in FeynRules [111,112], which generates UFO files [113] as output, which we then couple to MadGraph5. MadGraph5 outputs LHE event files, which will be further processed in different ways for ATLAS and the far detectors. These steps will thus be explained separately in the following subsections.
We note that the FeynRules model data base already contains two independent HNL implementations [114,115]. However, neither of these contains NROs. To include the operators of table I, we have therefore written our own model files, which are roughly based on the original HNL model implementation of Ref. [114].
The implementation including NROs is valid, strictly speaking, only for Dirac neutrinos since, as mentioned in its manual [93], MadGraph5 can not handle Majorana fermions in operators with more than two fermions. This problem can be circumvented by implementing some particular renormalizable SM extension, which includes the necessary scalar (or vector) fields, that will generate the operators under consideration in the infrared. As discussed in section II, the simplest model to generate the operators in table I is adding three different scalar LQs to the HNL extended SM.
We have written two of these model files, one for Dirac and one for Majorana HNLs. 10 We have checked that for Dirac neutrinos the EFT model file containing NROs and the HNL+LQ implementation give the same cross sections in the appropriate limits. This is shown in fig. 1, where we compare cross sections for HNL pair production calculated in EFT with those calculated in a model with a scalar leptoquark, S d . This example calculation chooses only the first generation couplings non-zero. We have repeated this check for other generation indices. As can be seen, the two calculations agree very well in the limit of large Λ = m LQ and small HNL masses. The differences between the two calculations can be sizeable, once m N approaches m LQ , as expected. 11 Since in the EFT total cross sections scale simply as Λ −4 , one can calculate LHE event files for large values of Λ (or large values of m LQ ) and scale the total event number to the desired value of Λ to be simulated. Note that this interchangeability is a necessary condition for us, in order to be able to apply our Majorana HNL+LQ implementation to a simulation of Majorana HNLs in EFT, where, by setting the leptoquark masses to 10 TeV in our numerical simulation, the HNL+LQ model is effectively reduced to the EFT.
In fig. 2, we compare cross sections for pair production of HNLs for Dirac and Majorana cases. The calculation is done with the HNL+LQ implementation and again with non-zero coupling to the first generation quarks only. For small HNL masses, the two calculations coincide. However, for m N larger than about 100 GeV numerically important differences show up. In order to understand this behaviour, we have calculated analytically the cross The results for Dirac and Majorana HNLs are given by We have also checked cross section calculations for different choices of operators. Figure 3 shows some examples. Since the LHC collides protons, one finds in general the cross section triggered by O QN to be larger than that induced by O uN , that in turn is larger than the one for O dN , as expected. In addition, we have calculated cross sections with all three N R pair operators switched on at the same time. Figure 3 shows an example with all couplings to N R set to one. This choice produces the largest cross section possible and thus represents the most optimistic assumption possible, as far as sensitivity is concerned. Typically, this choice leads to cross sections larger than the choice of only O dN , with c 11 dN = 1, by a factor of four or so. These two cases therefore span the "smallest" and the largest cross sections possible in our setup. 12

B. ATLAS
For the ATLAS experiment, our analysis will focus on the HNLs decaying to ejj. To enforce the HNL decays into ejj and to ensure numerical stability for very small decay width computation, we process the LHE file generated by MadGraph5 for the HNL decay in MadSpin [118]. The decayed LHE event files are then read by Pythia8 [119], which performs showering and hadronization. The displaced decay positions of N R into ejj can be extracted from the properly simulated input LHE files.
We generate 100 thousand events at each grid point in a two-dimensional scan in the m N -|V eN | 2 plane, for which we consider 56 values of m N from 5 GeV up to 6 TeV, and 175 values of |V eN | 2 from 10 −27 to 8 × 10 −3 , both in logarithmic steps. Event selection is then performed at the Pythia-event level, with a customized detector simulation to quantify the detector response to the physical objects of interest. Our displaced search design is inspired by the ATLAS inner tracker displaced-vertex (DV) search in Ref. [120], for the "DV + electron" channel. In this experimental search channel, photon triggers are used for final states involving displaced electrons. As these triggers have no cuts nor vetoes on tracks in the inner detector, displaced electrons can be found within these data sets. An electron candidate with p T > 120 GeV is required to have triggered the event in this channel [120]. 13 A similar search was proposed for displaced HNLs from Z decays in Ref. [57].
Our event selection starts by identifying electrons with |η| < 2.47 and p T > 120 GeV. A flat efficiency of 70% is assigned for electron identification. We then select events containing at least one reconstructed DV. For this reconstruction, we identify the N R displaced position at the truth-level and apply parameterized detection efficiencies as a function of DV invariant mass and number of tracks, as provided in the auxiliary material of the ATLAS 13 TeV search in Ref. [121]. 14 A reconstructed DV must lie within the acceptance of the ATLAS inner tracker. This means selecting displaced vertices with transverse position, r DV , to be within 4 mm < r DV < 300 mm and longitudinal position, z DV , to be |z DV | < 300 mm.
We require a displaced vertex to be made from at least four charged particle tracks. Each 13 As the ATLAS "DV + electron" analysis -as well as our displaced signal -does not possess a prompt isolated lepton, experimental triggers face the challenge to perform dedicated reconstruction techniques in order to enhance sensitivity to such LLP decay topologies. 14 For the purpose of our estimations in section IV, we assume these efficiencies remain the same at 14 TeV. track coming from the DV is required to have p T > 1 GeV and to be displaced. The latter is ensured by requiring |d 0 | > 2 mm, with d 0 = r · ∆φ defined as the approximate transverse impact parameter, with r corresponding to the transverse distance of the track from the interaction point (IP) and ∆φ being the azimuthal angle between the track and direction of the long-lived N R .
One of these displaced tracks must correspond to the electron that passed the trigger. As our electrons are reconstructed without an isolation requirement, we can match the truth-level index of the electron track with one of the displaced tracks.
The invariant mass of the DV, m DV , can then be reconstructed from the above selected tracks, and it assumes all the reconstructed tracks have the mass of a pion. The additional requirement of m DV ≥ 5 GeV is imposed in order to eliminate SM background from Bmesons. This last cut, together with the requirement of a large track multiplicity from the DV, defines a region where the signal is expected to be found free of backgrounds [120,121]. 15 In fig. 4, we show in two plots the individual and cumulative selection efficiencies of each cut for m N between 10 and 1000 GeV, for a Dirac HNL in a scenario with only the O 11 dN operator switched on. We have excluded the selection effect of the fiducial volume for these efficiencies, rendering the latter independent of the active-heavy mixing which has no impact on the other cuts than the DV fiducial volume one. We find a dip in the electron trigger efficiency at m N ∼ m W , where the transition between two-body and three-body decays into an on-shell and off-shell W -boson, respectively, takes place. The individual efficiencies of requirements on the number of tracks associated to the DV and on the DV invariant mass improve with increasing HNL mass, as expected. The parameterized DV efficiency cut in general passes most of the events for the shown mass range. As a result of these individual selections, the cumulative efficiencies tend to be enhanced for larger m N , except for the drop around the W -boson mass. Moreover, the m DV ≥ 5 GeV event requirement determines the lower mass reach.
The total number of signal events at the ATLAS is calculated with the following expression: where labels the efficiency of event selections including the ATLAS detector geometries, which depends on both the HNL mass and proper lifetime. Also, the cross section σ is a function of m N and the Wilson coefficients of the operators switched on, and Br(N R → ejj) depends on m N only. This efficiency includes detector acceptance and corresponds to the efficiency for reconstructing one displaced vertex in an event after all the selections (whereas the efficiencies shown in fig. 4 do not include the detector acceptance, i.e., equivalently, the fiducial volume requirement). Note that the factor of 2 arises from the fact that we simply require either of the two N R 's in each event to decay to ejj with a displacement while disregarding how the other HNL decays. Under the zero background assumption, 95% C.L. exclusion limits can be derived by requiring three signal events. These limits are provided and discussed in section IV. We note that constraints on mixing in the muon sector could also be obtained if the HNL decays to µjj. In this case, as the lepton coming from the DV is a muon, muon triggers can be used to efficiently record events, and we would expect similar exclusion reach for the muon case 16 (see for instance Ref. [57]). For mixing in the tau sector, we expect the sensitivity reach to be less powerful as a result of worse efficiencies and experimental difficulties in reconstructing (displaced) tau leptons.

C. Far detectors
Besides the default local experiments such as ATLAS and CMS, in recent years a series of far-detector (FD) experiments have been proposed with a distance of about 5-500 meters from the various IPs at the LHC: ANUBIS 17 , AL3X, CODEX-b, FASER and FASER2, MoEDAL-MAPP1 and MAPP2, and MATHUSLA. In particular, MoEDAL-MAPP1 and FASER have been approved. These experiments are supposed to be shielded from the associated IP by rock, lead, or other material, removing the SM background events effectively.
The other background sources include cosmic rays (especially for the on-the-ground experiment MATHUSLA), which can be eliminated by directional cuts. Consequently, in almost all the cases, essentially zero background can be assumed, and correspondingly the 95% C.L. exclusion limits can be derived by requiring three signal events. It is also worth mentioning that these experiments are proposed to be operated in different phases of the (HL-)LHC and IPs, and so the projected integrated luminosities vary. For instance, MAPP1 is to receive 30 fb −1 of integrated luminosity while FASER2 about 3 ab −1 .
The treatment of the sensitivity estimate at the far-detector experiments is somewhat different from that for the ATLAS experiment. In contrast to the ATLAS simulation, the decay of the HNLs is not simulated in this case. With MadGraph5 we generate 10 million events for 81 mass values from 0.1 GeV to 6 TeV in logarithmic steps. The LHE event files are fed into Pythia8 which simply includes the effect of initial state radiation and final state radiation. For 175 values of the mixing angle squared, |V eN | 2 , from 10 −27 to 8 × 10 −3 we then compute the decay probability of each simulated HNL in the fiducial volume (f.v.) of each far-detector experiment with the formulas below: where P [N R decay in f.v.] denotes the average decay probability of all the simulated HNLs (2 in each event) inside the fiducial volume of a far detector, k = 10 7 is the total number of the simulated events, and Br(N R → vis.) labels the HNL decay branching ratio into visible states, meaning all or some of the decay products are electrically charged (here we only consider the tri-neutrino final state to be invisible). P [N i R decay in f.v.] is the individual decay probability of the ith simulated N R , and it is computed with the help of the exponential decay distribution, kinematics (boost factor and traveling direction of the HNL), and the proper decay length cτ of the HNL. These expressions have to be calculated for each experiment with the corresponding cuts on the polar and azimuthal angles, as each of them is to be placed at a different location and thus has a different geometrical acceptance. Essentially P [N i R decay in f.v.] is calculated with the following formula: where L 1 is the distance from the IP to the point where the long-lived HNL reaches the detector, L 2 the distance the HNL travels through the internal space of the detector given its direction if it does not decay before leaving the detector chamber, and λ i is the decay length of the ith simulated HNL in the lab frame. P [N i R decay in f.v.] = 0, of course, if N i R would miss the detector volume. We emphasize that here, Eq. (20) is given only in order to showcase the essence of the computation of decay probabilities of the long-lived HNLs in a far detector. Of course, the far detectors considered in this work are proposed with various and even quite complicated geometries and relative positions from their respective IPs, for which one cannot directly apply Eq. (20). When we conduct the numerical simulation, we use more sophisticated and accurate formulas. We refer the reader to Ref. [53] and the references therein, for a summary of these far-detector experiments and the explicit expressions for the computation of P [N i R decay in f.v.].

IV. RESULTS
Based on the computation procedure as described in section III, we have estimated the 95% exclusion limits under zero background assumption (we will refer to them as "experimental sensitivities") on the effective operators containing two HNLs and two quarks O dN , O uN and O QN , for both ATLAS and future far detectors. For this purpose, we have treated the mixing with the active neutrinos |V αN | 2 with α = e, µ, τ , the HNL mass m N , and the operator coefficients c O /Λ 2 as independent parameters. For simplicity, we consider the case of n N = 1 generation of HNLs and focus on their mixing with the active neutrinos in the electron sector only, i.e., V αN = V eN .
In the theoretical framework considered in this paper, the HNL production at the LHC may occur either via the effective operators with a pair of HNLs or the mixing with the electron neutrinos (mediated by W and Z bosons). However, the decay of the HNLs can only proceed via the standard active-heavy mixing, since the HNL-pair EFT operators cannot induce the lightest HNL to decay. The HNLs can thus undergo either two-body or three-body (leptonic or semi-leptonic) decays, depending on their mass.
In fig. 5, we show the estimated experimental sensitivities at the 14 TeV LHC to a longlived Majorana HNL, in the plane |V eN | 2 versus m N . For our analysis we have considered different choices of operators. In the left panel, we focus on one single operator, O dN , choosing only the coupling with down quarks, c 11 dN , to be non-zero. In the right panel, we switch on all the three types of operators O dN , O uN and O QN , taking the same value for all the couplings of the first quark-generation combination, c 11 dN = c 11 uN = c 11 QN . These two choices lead to respectively the smallest and largest possible cross sections 18 for the HNL pair production at the LHC, and thus cover the most conservative and the most optimistic scenarios. For m N 5 GeV, the HNLs are produced dominantly by B-and D-mesons decays. These decays can proceed via mixing with active neutrinos as well as through d = 6 interactions. A detailed study of the interplay between the minimal scenario and the effective operators containing one N R has been performed in Ref. [53]. Inclusion of the N R pair operators in the meson decay rates goes beyond the scope of the present work and is to be performed elsewhere. Since we focus in our current numerical study on coefficients for quarks of the first generation only, we believe our figures will be affected by this approximation only in a minor way. Therefore, the inclusion of effective operators does not affect the LHC sensitivities to |V eN | 2 and in the region m N 5 GeV, we simply reproduce the results previously derived in the literature [49][50][51][52][53]. PS191 [24], JINR [25], and DELPHI [26], and the light gray band corresponds to the parameter space with the type-I seesaw relation assumed for active neutrino masses between 10 −3 and 10 −1 eV.
For m N 5 GeV, the HNL production is dominated by the effective operators. In this case, the mixing is only important for the HNL decays (but not for their production), and as a result, the sensitivity to the mixing grows substantially because of the much enhanced production rates of the HNLs. When the mixing and mass lead to a boosted decay length that falls well into the sensitive range of a detector, a large acceptance rate is expected. For larger m N , only a smaller mixing angle can roughly correspond to a fixed acceptance rate. As we can see in fig for an integrated luminosity of 3 ab −1 , for the most conservative (optimistic) scenario. Its lower mass reach is at about 5 GeV, which is due to the m DV > 5 GeV event selection we apply. ATLAS has a sensitivity on the mixing |V eN | 2 worse than MATHUSLA by almost two orders of magnitude, but it can reach larger masses of the HNL. This is mainly due to the distance at which MATHUSLA is planned to be located (∼ 100 m from the interaction point) which makes it sensitive to higher values of the lifetime, and in turn to smaller mixing |V eN | 2 and masses m N . As we can see, with the inclusion of the effective operators, these experiments can reach limits on the mixing |V eN | 2 several orders of magnitude better than the current experimental bounds, represented here with a dark gray area at the top of each plot in fig. 5. the experimental sensitivity reach to the values of mixing |V eN | 2 and mass m N required by the type-I seesaw mechanism to explain the neutrino masses. The future far detectors and ATLAS are complementary to each other in this regard. Figure 6 shows the estimated experimental sensitivities of both ATLAS and future far detectors in the |V eN | 2 versus m N plane, for one single effective operator O uN . We switch on here only the coupling with the up quark, c 11 uN . The left (right) plots are for a Dirac (Majorana) HNL 19 , for c 11 uN /Λ 2 = 1/(2 TeV) 2 . As expected in this case, the limits are slightly better than those for O dN (the upper left plot of fig. 5) and slightly worse than the limits shown in the upper right plot of fig. 5 (where all three operators for the first quark generation are switched on). The lower row of fig. 6 contains plots which are the zoomed-in version of those in the upper row, focusing on the large-mass regime. As we can see, the range of experimental sensitivity on the HNL mass m N is larger for the Dirac case than for the Majorana one. For instance, the mass reach of MATHUSLA, ANUBIS, and AL3X 19 For Dirac HNLs with mass below about 5 GeV, the limits are recast from the Majorana ones that exist in the literature [49][50][51][52][53], by multiplying |V eN | 2 with √ 2.
is at about 2.4 TeV for Dirac HNLs, compared to 2.0 TeV for Majorana HNLs. ATLAS, on the other hand, for an integrated luminosity of 3 ab −1 , can be sensitive up to masses of m N ∼ 3.0 (3.6) TeV for a Majorana (Dirac) HNL. This difference in the mass reach between Dirac and Majorana HNL scenarios is due to the different HNL pair production cross sections in the two cases. In particular, this difference is more perceptible for larger m N (cf. eqs. (15) and (16), fig. 2 and the related discussion in section III A). Overall, the shape of the main results shown in figures 5 and 6 can be qualitatively understood as the final efficiency being bounded by the cases when the HNLs are decaying either too promptly (the upper right corner) or too far away (the lower left corner), both outside the detectors' acceptance. In addition, the upper mass reach is limited by the production cross sections.
We would like to note that Big Bang Nucleosynthesis can place further constraints on the active-heavy mixing angle in the low mass region m N 2 GeV (see, e.g., Refs. [123][124][125]). However, these constraints rely on the evolution of HNLs in the early Universe, discussion of which goes beyond the scope of our work. Thus, we remain agnostic about them and do not show the corresponding excluded regions in our figures.
Finally, fig. 7  Currently the LHC could potentially reach values of Λ ≈ 3 TeV when reinterpreting 13 TeV prompt searches. For a reinterpretation of the prompt same-sign dilepton search performed by CMS in Ref. [27], see appendix B.

V. SUMMARY
In this work, we have studied the phenomenology of heavy neutral leptons (HNLs) as long-lived particles (LLPs) in the framework of an effective field theory (EFT). Concretely, we have studied the effects of d = 6 four-fermion operators with a pair of HNLs and a pair of quarks in the EFT of a SM extension with HNLs as "right-handed neutrinos", N R SMEFT. We considered three different types of these operators and have studied scenarios with only one or all of them being present at the same time. While quantitatively the results depend on the assumptions made about type of operator and/or Wilson coefficient present, qualitatively all HNL pair operators behave similarly. Operators with two HNLs and two quarks lead to production cross sections at the LHC which are not suppressed by the small mixing of the HNLs with the active neutrinos. Instead, cross sections are proportional to Λ −4 , where Λ is the energy scale of the non-renormalizable d = 6 operators. HNL pair operators will not cause (the lightest) HNL to decay. Thus, HNLs decay only via their mixing with the active neutrinos, V 2 . This scenario leads to one very important change in phenomenology with respect to previous works, that considered HNLs at the LHC which were produced via charged (and neutral) current diagrams induced by the mixing of the HNLs only. Namely, the total signal event number for the different experiments in our setup scales at the smallest V 2 that one can probe only as V 2 , instead of the usual V 4 .
We have estimated the sensitivity range of various LHC experiments in this setup: ATLAS and a series of proposed "far detectors", AL3X, ANUBIS, CODEX-b, FASER, MATHUSLA, and MoEDAL-MAPP. Our main result is that for Λ lower than roughly Λ 10 − 15 TeV, depending on the operator, much larger HNL masses and much smaller mixing angles, V 2 , can be probed than in the "standard" case, where both, production cross section and decay lengths, are determined by the mixing angle (and the HNL mass) only.
We have also briefly discussed some ultraviolet completions for the considered d = 6 operators: leptoquarks (LQs) and a possible Z . While we presented our results only in terms of EFT parameters, it is clear from the numbers found in our simulation that the reach in LQ (or Z ) mass in these "indirect" searches will be much larger than for direct on-shell production at the LHC.
We have also highlighted differences in the production cross sections for Dirac and Majorana HNLs and shown both, analytically and numerically, that cross sections at HNL masses larger than roughly 100 GeV are different in the two cases. The relative suppression of Majorana HNL cross section with respect to the Dirac one can not, however, be easily used to distinguish the two cases: velocity distributions for Majorana and Dirac HNLs are different at large HNL masses, but measuring these would require a rather large event number. However, this change in cross section results in a smaller sensitivity reach of the different experiments for Majorana HNLs at the largest HNL masses that can be probed at the LHC. Rather, to distinguish between Majorana and Dirac HNLs one should search for events in ATLAS, in which both N R 's decay inside the detector. If such an event can be identified for Majorana HNLs the ratio of same-sign to opposite-sign dilepton events should equal 1. This could be checked with just very few events, in principle.
It is possible to require reconstruction of both DV's in each signal events in ATLAS (not in the far detectors, however). This would also be more effective at suppressing background and hence allow for looser event selections cuts, enhancing the event selection efficiencies. However, it would also reduce the signal event number, since the probability of a N R decaying inside the pixel detector is always smaller than one and in the extremes of parameter space that one can probe actually much smaller than one. As a result, while observing double DV events would be very interesting for establishing the Dirac/Majorana nature of the HNLs, searching for such events will not lead to the possibility to explore new parts of the parameter space of N R SMEFT.
Finally, we close by mentioning that there are also d = 6 operators with only one HNL. Parts of the parameter space where such operators may be probed in the far detectors have already been studied in Ref. [53].

Appendix A: A timing cut at ANUBIS
In this appendix, we discuss briefly the effect of a potential timing cut at the ANUBIS experiment on the expected exclusion limits for our setup. ANUBIS [40] has been proposed to be installed inside one of the service shafts above the ATLAS detector. The relative proximity of ANUBIS to ATLAS would allow ANUBIS to trigger on the readout of ATLAS. Given the close distance from the IP, however, it is not as straightforward to implement background vetos as in some of the other proposed far-detector experiments; some discussion about backgrounds is given in the original proposal [40].
One of the possible ways to reduce background events is to make use of timing information. 20 Since the hard collision can be timestamped, if a DV is formed significantly later than expected, the event can be rejected with this strategy, such that background events from e.g. beam-gas and beam-collimator interactions could be reduced to a large extent. Such a timing cut is based on the idea that light LLPs (masses of a few GeV) are produced at the LHC with velocities β 1, where c labels the speed of light. In the EFT setup, discussed in this paper, however, HNLs with much larger masses can be probed and their velocity distributions will depart from β = 1, leading to a significant delay in arrival times of the events in ANUBIS.
Here, for a rough estimate of the size of the effect, we consider two possible acceptance time windows: 1 ns and 2 ns. Note that these effectively correspond to LLP speeds of 0.99 · c and 0.98 · c for a distance of about 30 m which is roughly the distance of ANUBIS from the ATLAS IP. As a benchmark, we take Dirac HNLs with the single NRO O 11 dN . Results for the other operators will be qualitatively very similar. Figure 8 shows two plots to exemplify the effect of a timing cut. The left plot shows various curves for the ratios of the numbers of the HNLs traveling inside the solid angle coverage of the bottom of ANUBIS with the speed above certain thresholds to the total number of the produced HNLs, as a function of m N . We find that for m N 200 GeV the fraction of HNLs in the direction of ANUBIS with a speed larger than 0.99 · c drops below 10% of the total number already, and for a mass close to 1 TeV the ratio even drops down to the level of 10 −4 . The effect of this loss of events on the sensitivity reach is demonstrated in the right panel. For m N above about 200 GeV, a timing cut will reduce the exclusion power by a significant amount. For instance, if we require 3 signal events, excluding HNLs with a speed smaller than 0.99 · c, the reach in m N shrinks from about 2.1 TeV to only 360 GeV and in |V eN | 2 from 5 × 10 −23 to 3 × 10 −20 .
The plot on the right also shows curves for two different assumptions of excluded events, 3 and 30 signal events, respectively. We conclude that for the largest HNL masses it could  be advantageous to use a looser timing cut, since larger parts of parameter space could be probed this way, despite larger backgrounds. We close this discussion by noting that a more detailed MC simulation of backgrounds in ANUBIS would be necessary in order to determine the optimal cuts, which is beyond the scope of our work. thus upper limits on V 2 eN and V 2 µN are presented as a function of HNL mass, m N . 21 Note that this search explicitly requires same-sign leptons and thus is valid only for Majorana HNLs. One can expect limits on Dirac HNLs to be significantly weaker in such a search, since opposite-sign leptons have much larger SM backgrounds.
The search [27] assumes the HNLs decay promptly. A lower limit on Λ, calculated from this data, will therefore be valid only if V 2 αN is large enough, such that a sufficient number of events decay within the cuts used by CMS. We will quantify the range of V 2 αN as function of m N , for which this is the case, in more details below. Let us first concentrate on how to estimate limits on Λ.
Since CMS decided to present their results only in the plane (V 2 αN , m N ), the first step in our reinterpretation is to calculate the corresponding number of excluded signal events as function of m N . In their modeling, CMS [27] takes into account two types of Feynman diagrams for the production of HNLs: (1) s-channel Drell-Yan (DY) process, mediated by a W -boson and (2) "photon-initiated" process, which CMS also calls VBF (vector-bosonfusion) channel. The importance of the latter at large HNL masses was pointed out in Refs. [128,129]. This is due to the fact that VBF diagrams contain a W -boson in tchannel, which has a softer m N dependence than the Drell-Yan diagram at proton colliders. 22 VBF diagrams do depend very sensitively on the photon content of the proton. In a recent paper, Manohar et al. [130] presented a set of PDFs that improve the calculation of the photon density over previously existing PDFs. CMS uses in the calculation the LUXQED17 PLUS PDF4LHC15 NNLO 100 PDF set [130], which we also use for our reinterpretation for consistency.
CMS provides sample values of cross section for DY and VBF for three choices of m N = 40, 100 and 1000 GeV. Using MadGraph5 [91][92][93], we recalculated the cross sections for these points as a cross-check of the reliability of our conversion procedure. While for m N = 40 GeV we reproduce very well the numbers quoted by CMS, for larger m N our results depart from the numbers given by CMS. At m N = 1 TeV, our VBF (DY) result is ∼ 75% (80%) larger than those by CMS. We can only speculate that CMS uses some additional cuts in the run card.dat of MadGraph5, that we were unable to locate in the CMS paper. We decided therefore to use the difference of our cross section results from those of CMS as an estimate of the uncertainty for the total event number extraction in our reinterpretation procedure. 21 Also the ATLAS collaboration has published similar searches. However, [126] is based on √ s = 8 TeV data, and thus no longer competitive. The analysis of [127], on the other hand, concentrates on a leftright extension of the SM in its analysis. The changes in cross section and kinematics for on-shell W R relative to the expectations for N R SMEFT, however, would make a reinterpretation of this paper even more unreliable than our estimates for [27]. We therefore decided to concentrate on the CMS search here. 22 The final state of the VBF diagrams contains an additional jet, which will be mostly central. The same final state can also be generated by "gluon-initiated" diagrams. These are numerically more important than the "photon-initiated" diagrams, except for m N in excess of roughly 1 TeV. (The exact number depends very strongly on the PDF.) However, CMS does not take into account these diagrams in their modeling, so we have to exclude them in our reinterpretation too. The upper limits on V 2 eN plus the calculated cross sections allow us to estimate the number of excluded signal events (as a function of m N ). This number, however, includes the detector efficiencies indirectly. So, in order to extract limits on, for example, Λ in our EFT setup we have to assume that detection efficiencies are the same for both models. While this is probably a good approximation for small m N , say below the W −threshold, our operators produce pairs of HNLs, which at large m N affects the kinematics, and thus efficiencies. Given the uncertainty in total cross sections discussed above, we believe this to be currently a less important source of error for our estimate and will simply assume the efficiencies are the same in both cases.
We then show in fig. 9 to the left limits on Λ for the operator O dN as a function of m N . As in the main text, we consider only the first-generation quarks throughout this appendix. The red line is our estimate for the sensitivity of Ref. [27]. Since this search is background dominated, we can estimate the future sensitivity by scaling with the square root of the luminosity increase. This results in the blue line, representing the final high-luminosity LHC result. The corresponding coloured bands are our estimate of the uncertainty of our conversion procedure, based on the cross section uncertainties. The bands are made symmetric, assuming for simplicity that our cross section has a symmetric error of the order of the difference of our calculation relative to the CMS numbers.
As mentioned above, the search by CMS assumes that the HNL decays promptly. The cuts on the impact parameter used by CMS correspond to 0.1 (0.4) mm in the transverse (longitudinal) coordinates. We therefore simply require HNLs to decay with lengths smaller than L exp = 0.1 mm. This allows us to estimate the range of V αN versus m N that could be probed, once we fix Λ (and the type of operator). An example for O dN and Λ = 2 TeV is shown in fig. 9 on the right. Because of the fast decays for HNLs with masses larger than about 100 GeV, very small values of V αN could be probed this way, if Λ is small enough to yield testable production cross sections. We stress, however, that already for Λ in excess of 3 TeV, currently no new constraints on V αN could be set, beyond those already found in the CMS analysis. The results shown in fig. 9 are valid for α = e. Results for α = µ should be quite similar. Also, in the example calculation we used O dN for production of HNLs. Results for the other HNL pair operators are qualitatively similar.
We close this discussion by stressing again that many HNL searches at the LHC could be used to place limits on N R SMEFT parameters. Our reinterpreation of the CMS search [27], however, carries a large uncertainty and reinterpretation of other searches will be similarly limited. We believe only the experimental collaborations will be able to make reliable searches for our EFT scenario. The main issues we found in this reinterpretation are: (i) our conversion of the limits presented in Ref. [27] into limits on cross section times branching ratio has a surprisingly large error. We therefore urge our experimental colleagues to present the "model-independent" cross section times branching ratio limits in addition to the model-dependent final results in future publications. And, (ii) for an improved estimate a better understanding of the differences in efficiencies for HNL detection in different models would be needed. However, this will require a complete Monte Carlo simulation, which is beyond the scope of our present paper.