Long-lived heavy neutral leptons at the LHC: four-fermion single-NR operators

Interest in searches for heavy neutral leptons (HNLs) at the LHC has increased considerably in the past few years. In the minimal scenario, HNLs are produced and decay via their mixing with active neutrinos in the Standard Model (SM) spectrum. However, many SM extensions with HNLs have been discussed in the literature, which sometimes change expectations for LHC sensitivities drastically. In the NRSMEFT, one extends the SM effective field theory with operators including SM singlet fermions, which allows to study HNL phenomenology in a “model independent” way. In this paper, we study the sensitivity of ATLAS to HNLs in the NRSMEFT for four-fermion operators with a single HNL. These operators might dominate both production and decay of HNLs, and we find that new physics scales in excess of 20 TeV could be probed at the high-luminosity LHC.


Introduction
Interest in long-lived particles (LLPs) has grown largely in the last few years [1][2][3]. Many models for LLPs have been discussed in the literature, most of which are motivated by either dark matter or neutrino masses. Heavy neutral leptons (HNLs) are the prime example for LLPs connected with the neutrino masses. HNLs are Standard Model (SM) singlet fermions that couple to SM particles via their mixing with active neutrinos.
The minimal model that can realize this effective setup is the seesaw mechanism, in which right-handed Majorana neutrinos, N R , are added to the SM particle content [4][5][6][7][8]. However, many SM extensions, that aim to explain observed neutrino data [9] (see also refs. [10,11]) via electroweak scale variants of the classical seesaw [12][13][14][15], do not include only HNLs. For example, right-handed neutrinos appear necessarily in left-right (LR) symmetric extension of the SM as the neutral component of the right-lepton doublet [16,17]. If the additional non-SM states, such as the W R and Z in the LR model, have masses which are too large to be produced on-shell at the LHC, their effects on HNL phenomenology is best treated in effective field theory (EFT).

JHEP01(2022)044
Phenomenological interest in this EFT is motivated by the future upgrades of the LHC on one side and the improvement in the sensitivities of low-energy experiments on the other.
Effective interactions of d ≤ 6 are the most interesting from a phenomenological point of view. There are two d = 5 operators involving N R . Their phenomenology has been studied in detail in refs. [22,40,41]. The d = 6 operators containing N R can be divided into two classes: (i) operators with two fermions and bosons and (ii) four-fermion operators. The second class, in turn, can be partitioned into operators with two N R 's and operators with a single N R . 2 The LLP phenomenology of pair operators has recently been studied in ref. [39]. Here, we will concentrate on operators with a single N R . The phenomenology of single-N R operators is decidedly different from that of pair operators. First, pair operators do not by themselves lead to decays of (the lightest) N R . Instead, for these operators N R decays are controlled by the mixing with active neutrinos. This is different from the single-N R operators, which will usually dominate the decay length of the HNLs in those parts of parameter space where the operators are large enough to dominate N R production. Thus, the parameter space that can be explored for these two types of operators is very different, see section 4. Second, pair operators do not produce prompt charged leptons, except in the parameter region where the decay length of the N R is so short that the lepton from a N R decay is confused with a charged lepton produced directly from pp collisions at the interaction point (IP). In all lepton-number-conserving single-N R operators, on the other hand, N R 's are accompanied by a prompt lepton (either a neutrino or a charged lepton). This affects the search strategy for the different operators.
Ref. [38] studied single-N R operators for various proposed LLP "far" detectors, such as MATHUSLA [3,42,43], CODEXb [44], AL3X [45], FASER [46], and ANUBIS [47], as well as ATLAS, for HNLs produced from charm and bottom meson decays and hence with mass below 5 GeV. 3 In addition, ref. [51] very recently worked on phenomenology of the same set of single-N R operators associated with the third-generation leptons at Belle II, for HNLs produced from τ lepton decays. For these reasons, in our numerical simulation we concentrate on ATLAS for heavier HNLs, and a short discussion will also be given for the expectations for CMS (see section 4).
The rest of this paper is organized as follows. In the next section, we will discuss briefly N R SMEFT at d = 6. This section also entails a short discussion on how the single-N R operators could be the low-energy remnant of some leptoquark or two Higgs doublet models. Section 3 discusses the details of the simulation we perform for the ATLAS detector. In section 4, we present our numerical results. First, we discuss again briefly the minimal case, in which HNLs are produced and decay via mixing only. While this was previously done by some of us in ref. [52], we now also simulate the expectations for HNLs coupled to τ 's, including both neutral and charged currents leading to more realistic estimates for the future ATLAS sensitivities. We then present our results for the different single-N R operators. Cross sections and decay lengths depend on both, operator type and generation 2 There is also a lepton-number-violating operator with four NR's, but it requires at least two generations of HNLs. 3 For the expectations for these experiments in the minimal HNL scenario with only active-sterile neutrino mixing, see for example refs. [48][49][50].

JHEP01(2022)044
indices in the SM sector. For the first generation of SM quarks, sensitivities will reach new physics scales in excess of 20 TeV at the high-luminosity LHC. We then close with a short summary of our results.

Effective interactions
In this section, we briefly introduce the N R SMEFT, focusing on the operators of interest for the current work. If HNLs with masses below or around the electroweak scale exist in nature, the effects of new multi-TeV physics at much smaller energies can be systematically described in terms of an EFT built out of the SM fields and N R . At renormalizable level, in addition to the SM operators, there are a Majorana mass term for N R and a d = 4 operator describing the fermion portal: are the Wilson coefficients, and the second sum goes over all independent interactions at a given dimension d. At d = 5, in addition to the renowned Weinberg operator composed of L and H [53], one finds two more operators that involve N R [21,22].
At d = 6, in addition to the pure SMEFT operators [54], there are five operators involving two fermions (at least one of which is N R ) and bosons, eleven baryon and leptonnumber-conserving four-fermion interactions, one lepton-number-violating operator, and two operators that violate both baryon and lepton number [24]. 4 In the present work, we are interested in the effects of the lepton-number-conserving four-fermion interactions containing one N R and three SM fermions. We list them in table 1. The effects of the fourfermion operators containing a pair of HNLs and a pair of quarks have been investigated in detail in ref. [39].
The single-N R operators including quarks can lead to enhanced HNL production cross section at the LHC, but they also trigger the decay of N R to a lepton and two quarks. The total decay width of the N R 's depends on the operator. Neglecting the masses of the lepton and light quarks, the partial decay width to charged leptons plus quarks is given by 3)

JHEP01(2022)044
Name Table 1. Lepton-number-conserving four-fermion single-N R operators. For each operator structure, we provide the number of independent real parameters for n N = 1 and n N = 3 generations of N R . The operator in the last row is purely leptonic, and thus, it does not contribute to the HNL production at the LHC. 3) applies to Dirac neutrinos. For Majorana neutrinos, one has to add also the charged conjugated channels, leading to another factor of 2 for the widths.

Ultra-violet completions for four-fermion single-N R operators
The single-N R operators of interest can be generated in ultra-violet (UV) complete models containing heavy scalars or vectors. Here, we do not aim to provide a complete classification of such UV completions, but rather give a few examples. In what follows, we consider scalar leptoquarks and an inert SU(2) L doublet scalar. A catalog of models with scalar and vector leptoquarks generating four-fermion operators involving one or two N R 's and quarks can be found in ref. [30]. The operator O duN e can arise from a model with a scalar leptoquark S d having the gauge quantum numbers of the down quark, cf. table 2. The interaction Lagrangian of S d is given by Upon integrating out S d , the operator O duN e is generated with the tree-level matching condition for the Wilson coefficient c duN e given in the last column of table 2. 5 Analogously, a scalar leptoquark S Q with the quantum numbers of the SU(2) L quark doublet can lead We note that the first terms in eqs. (2.4) and (2.5) also generate the N R pair operators where q = d R and q = Q, respectively, cf. ref. [39].

JHEP01(2022)044
Heavy scalar   [58,59], we have checked that both cases lead to the same single-N R production cross section. We note that for N R pair production triggered by the four-fermion operators with two N R 's, the cross section is different for Dirac and Majorana HNLs, especially for values of m N 100 GeV at LHC energies (we refer the interested reader to section 3.1 of ref. [39]). The fact that the HNL nature does not affect the production triggered by the four-fermion single-N R operators allows us to implement these operators directly in FeynRules for Dirac HNLs and use the resulting UFO model file in MadGraph5. (Recall that MadGraph5 can not handle Majorana fermions in operators with more than two fermions, cf. section 3.1 of ref. [39].)

Simulation details
Our signal topology contains a prompt lepton and a displaced vertex (DV) stemming from the N R decay to leptons and quarks. Our stage to reconstruct such a signature is the ATLAS detector, specifically its inner tracker, as it has the capability to reconstruct vertices displaced from the IP by few millimeters to tens of centimeters. Our analysis strategy builds up on an earlier work [52] and is inspired from ATLAS multi-track displaced searches [60,61].
We consider the collision process pp → N l with l = e, µ, τ , at √ s = 14 TeV at the highluminosity LHC with an integrated luminosity of 3 ab −1 . We generate LHE events with displaced information at the parton level with MadGraph5, which are read by Pythia8 [62] JHEP01 (2022)044 for showering and hadronization. Our detector simulation is based on a custom made code within Pythia8, where we first reconstruct isolated prompt electrons, muons, and taus (with help from FastJet [63]), taking into account detector acceptance, resolution, and smearing on their transverse momenta (for details, see ref. [52]). After selecting events with a prompt lepton, the displaced vertex reconstruction starts by selecting tracks 6 with p T > 1 GeV and a large impact parameter, d 0 , defined as d 0 = r trk × ∆φ. Here ∆φ corresponds to the azimuthal angle between the track and the direction of the long-lived N R , and r trk corresponds to the transverse distance of the track from the origin. We require |d 0 | > 2 mm.
As we have access in simulation to truth-level Monte Carlo information, we also identify the truth N R decay positions in the transverse and longitudinal planes, namely, r DV and z DV , respectively. An additional step (with respect to ref. [52]) of the vertex reconstruction implemented in this work is the requirement that r trk − r DV < 4 mm. It is not always the case that the "starting" point of the displaced track matches the displaced vertex position. This is more evident in the case where we have a tau produced from the N R displaced decay, as taus also have an additional displacement. 7 With this requirement, we emulate what an experimental displaced-vertex reconstruction would do when fitting nearby displaced tracks to a common origin [60]. This will lead to an additional reduction in efficiency when reconstructing displaced vertices containing taus. Nevertheless, it is a more realistic (and optimistic) approach than what was done in ref. [52] to handle heavy neutrino decays to taus (see section 4 below).
After selecting optimal displaced tracks, we demand displaced vertices within the AT-LAS inner tracker acceptance, namely, 4 mm < r DV < 300 mm and |z DV | < 300 mm. Further cuts are applied on the number of high-quality tracks coming from the DV, N trk , and its invariant mass, m DV , assuming all tracks have the pion mass. More concretely, we require N trk > 3 and m DV ≥ 5 GeV. As detailed in refs. [60,61,64], these last two cuts ensure that we are in a region where signal is expected to be found free of backgrounds including B-mesons. Further detector response to DVs is quantified by applying the 13 TeV ATLAS parameterized efficiencies [61] as a function of DV invariant mass and number of tracks, where we assume these will remain the same at 14 TeV.

Numerical results
Based on the computational procedure described in the previous section, we have estimated the experimental sensitivities (95 % confidence level (C.L.) exclusion limits under the assumption of zero background) of searches for long-lived HNLs at the ATLAS detector for two different theoretical scenarios. The first is the minimal scenario in which only right-handed neutrinos, N R , are added to the particle content of the SM and renormalizable interactions are assumed. In this case, the HNLs interact with the SM particles 6 A track in our simulation is a final state charged particle. These come from the decays of NR and can correspond to an electron, a muon, or a charged particle coming from the hadronization of quarks or from tau decays. 7 The proper decay distance of tau leptons is cτ = 87.1 µm. This will lead, for example, to decay distances of γcτ ∼ 5mm at 100 GeV.

JHEP01(2022)044
only through the mixing with the active neutrinos, V lN , with l = e, µ, τ . In the second theoretical scenario we consider N R SMEFT containing non-renormalizable interactions of N R with the SM. In this case, both production and decay of the HNLs can be mediated by the single-N R effective operators under consideration. In all of the plots below, we assume a 3 + 1 scenario where the HNL mixes dominantly with only one active neutrino flavor at a time. We assume only one HNL is kinematically relevant.

Minimal scenario
In the minimal scenario, the relevant parameters are the mass of the HNL, m N , and the mixing of the HNL with the active neutrinos, V lN , which we have treated as independent parameters. The HNLs are produced from the decays of on-shell W -bosons into a lepton and an HNL associated with a charged lepton, pp → W → lN , via the HNL mixing with the active neutrinos. The decay of the HNLs occurs also via the mixing with the active neutrinos, through both charged and neutral SM currents, N → l(ν)jj. For the minimal scenario we use the FeynRules implementation for HNLs of ref. [65]. Figure 1 shows the region, in the plane |V lN | 2 vs. m N , where a displaced-vertex search at the ATLAS detector for the center-of-mass energy 14 TeV, and with the selection criteria discussed in section 3, may have sensitivity to the minimal scenario. As can be seen in this figure, the sensitivities in |V eN | 2 and |V µN | 2 are rather similar and can reach values down to |V lN | 2 ∼ 10 −9 for m N ∼ 30 GeV, with 3 ab −1 of integrated luminosity. On the other hand, in the case of mixing with the tau neutrinos, ATLAS can reach values of the mixing parameter down to |V τ N | 2 ∼ 5 × 10 −9 for m N ∼ 20 GeV with 3 ab −1 . Figure 1 compares our limits with the current experimental bounds for this model, represented by the dark gray area at the top of each plot. These constraints were obtained at the following experiments: ATLAS [66], CMS [67,68], DELPHI [69], and LHCb [70,71]. As we can see, our forecast limits can reach values of the mixing |V lN | 2 several orders of magnitude smaller than current experimental bounds.
It might be interesting to compare the forecast limits to theory expectations. In seesaw type-I, one naively expects |V lN | 2 m ν /m N (10 −12 − 10 −11 ) for values of m N in the range we are considering in this work. Larger values of |V lN | 2 are possible allowing finetuning in parameters. A more natural model for having |V lN | 2 in the range accessible for ATLAS/CMS might be the inverse seesaw [12]. In this variant of the seesaw, |V lN | 2 is given by |V lN | 2 m ν /µ, with µ being the lepton-number-violating parameter of the inverse seesaw model. With µ supposedly being a small parameter, when compared to m N , mixings in this model can be easily as large as the experimental limits.
As mentioned above, the same search strategy for long-lived HNLs was previously proposed by some of us in ref. [52]. One of the differences between our current and previous calculations is the center-of-mass energy at the LHC, which now is taken as 14 TeV (previously in ref. [52] we used 13 TeV). Perhaps more important is the fact that in the present paper, our numerical calculations used more statistics, which allowed us to obtain much smoother contours for our limits, which led to a slight increase in the ranges shown. Moreover, in the case of mixing with taus, our current limits are more sensitive than the previous ones calculated in ref. [52]. The reason for this difference is that in ref. [52], JHEP01(2022)044  [67,68], DELPHI [69], and LHCb [70,71].
we only considered neutral currents in the decay of HNLs that coupled to taus (i.e. we ignored a tau lepton coming from the displaced vertex), whereas now, we have included both charged and neutral currents in our calculations, making our limits more realistic for the case of the mixing with the tau neutrinos and comparable with the sensitivity reach projected with other proposed strategies (see for instance ref. [72]).

Four-fermion single-N R operators
In the second theoretical scenario, we consider the four-fermion single-N R operators in the N R SMEFT. We estimate the experimental sensitivity of our displaced search to a long-lived HNL at the ATLAS detector. Here, we take the coefficients of the operators c O /Λ 2 and the mass of the HNL, m N , as independent parameters. In this scenario, both the production and the decay of the HNL can be dominated by the same operator O, unlike the case of effective operators with two HNLs [39], where the pair-N R operators dominantly induce the HNL production, but the decay of the HNL still proceeds only via mixing with the active neutrinos. For the EFT scenario, in our analysis we have assumed that the contributions to the production and decay of the HNL from its mixing with active neutrinos V lN are subdominant and negligible compared to the effective operators' contributions. For mixing angles smaller than |V lN | 2 10 −9 , this assumption is always fulfilled.
The production of the HNLs considered in our analysis, pp → lN , is always accompanied by a prompt charged lepton -an electron, muon, or tau, depending on the flavor JHEP01(2022)044 structure of the effective operator considered. The presence of this charged lepton is important in our analysis as it is used to trigger the signal, as discussed in section 3. The decay of the HNLs will occur via the same operator leading to two jets and one neutral or charged lepton, N → l(ν)jj. The production cross sections of the HNLs will depend on the type of quarks that the respective operator includes. In our analysis, we have only considered effective operators with quarks of the first two generations.
In figure 2, we show the experimental sensitivity of the ATLAS detector to a long-lived HNL in the Λ vs. m N plane. In our analysis, we have considered the contributions of one operator at a time, setting the value of the corresponding operator coefficient c O = 1, and the rest of the operator coefficients to zero. In figure 2, we have considered only operators with quarks of the first generation. Note that the numbers in the superscript of e.g. c 1112 duN e refer to the first-generation quarks (d and u), the lightest N R and the second-generation charged lepton (the muon). As can be seen in this figure, for an integrated luminosity of 3 ab −1 , ATLAS can reach values of the new physics scale up to (and above) Λ ∼ 20 TeV for masses m N 50 GeV in the case of operators with an electron or muon. In the case of operators with a tau lepton, ATLAS can reach Λ 10 TeV at masses m N of 10's GeV. It is worth mentioning that our limits start at m N 5 GeV. The reason is the kinematic cut at m DV ≥ 5 GeV imposed in the selection criteria. This cut is necessary to remove the SM background coming from B-mesons, as discussed in section 3. We also note that the JHEP01(2022)044   , d) and (c, s)). We therefore do not show results for these cases explicitly. Operators with third-generation quarks have not been considered in this work, since they will require special treatment (i.e. tagging).
We also note that figures 2 and 3 have been calculated for Dirac HNLs. As mentioned above, production cross sections for single-N R operators are the same for Dirac and Majorana HNLs, while the half-lives for Majorana HNLs are smaller by a factor of two. The sensitivity regions for Majorana HNLs therefore differ slightly from the regions shown in the figures. We do not repeat the plots for the Majorana case and instead opt for a short explanation of the differences. First, the maximal value of the HNL mass, to which this kind of search is sensitive is determined by the smallest decay length that is accessible in JHEP01(2022)044 the experiment. Since the decay width scales as m 5 N , for a Majorana HNL the largest HNL mass accessible is a factor (1/2) 1/5 0.87 smaller than in the Dirac case. Second, the maximal value of Λ reached in our sensitivity curves is essentially determined by the total cross section (times luminosity). Since cross sections are the same for Dirac and Majorana HNLs, this maximal value of Λ does not change for Majorana HNLs. Finally, in the regime where the decay lengths are large (i.e. for large values of Λ at values of m N smaller than the one where the maximal value of Λ is reached), the event number depends linearly on the half-life, while both the cross section and the decay width scale as Λ −4 . For Majorana HNLs, in this part of the parameter space, slightly larger values of Λ are accessible than for the Dirac case, i.e. an increase by roughly a factor 2 1/8 1.09.
Let us briefly comment on existing limits. First of all, our choice of switching "on" always only one Wilson coefficient at a time guarantees that there is no new source of lepton (or quark) flavour violation. It is well known that ultra-violet completions, such as the leptoquark models we discuss in section 2.2, are strongly constrained by searches for lepton-flavor-violating processes. These will put lower bounds on m LQ Λ that are much stronger than anything achievable in direct searches [73,74]. Thus, in all accelerator searches it is customary to assume that new resonances, such as leptoquarks, couple only to one SM fermion generation at a time. Direct searches for leptoquarks from pair and single leptoquark production have been performed by both, CMS [75][76][77] and ATLAS [78,79]. The best limits approach now m LQ 2 TeV. Thus, the long-lived particle search discussed in this paper, will probe so-far uncharted parts of parameter space. Let us also mention, that since d = 6 operators have Wilson coefficients of the form c/Λ 2 , limits on Λ will scale proportional to √ c. Thus for c 10 −2 our search will no longer probe values of Λ not already excluded by direct LHC searches.
We will close this discussion with one additional comment. Our simulated analysis focused on the ATLAS detector and its reconstruction capabilities to displaced vertices inside the inner tracker, starting from 4 mm in multi-track searches [60,61]. Relaxing this requirement to decay distances below 4 mm (both in d 0 , r DV and z DV ) will allow to extend the reach in parameter space towards larger HNL masses. Of course, with the loosening of these cuts we may depart from the zero background case assumption, and a detailed study on the multi-track search backgrounds would be needed, which goes beyond the scope of the present work. Nevertheless, past displaced lepton searches -whose tracks are fitted to a common vertex -at CMS [80] could probe transverse decay lengths starting from ≈ 200 µm. 8 In addition, a recent 13 TeV CMS search [81] demonstrates that lepton tracks with |d 0 | > 0.1 mm are displaced enough to be considered for analysis. Finally, the recent CMS note on an HNL search with an explicit displaced vertex requirement does not even demand a constraint on the DV minimal distance [68]. This provides feasibility to experimentally go below the 4 mm threshold. We stress that an improvement of the displaced vertex search towards smaller decay lengths by such a larger factor (up to 40 for 0.1mm) would allow to test HNL masses larger by a factor 2 w.r.t. the values in our JHEP01(2022)044 figures, i.e. extend the searches from m N 50 GeV to roughly 100 GeV. We hope that this large potential gain motivates the experimental collaborations to study the lowering of the transverse cuts in displaced vertex searches to the sub-millimeter range.

Summary
The Standard Model (SM) effective field theory (EFT) extended with sterile neutrinos, also known as the N R SMEFT, provides a framework to systematically study sterile neutrinos associated with a high new-physics (NP) scale in ultra-violet complete models beyond the SM. In the N R SMEFT, high-scale NP effects are encoded in the so-called Wilson coefficients of non-renormalizable operators at different mass dimensions. Higher-dimensional operators involving N R can have either one, two, or four sterile neutrinos, and may conserve or violate lepton number, or else both lepton and baryon numbers.
In this work, we have focused on lepton-number-conserving four-fermion single-N R operators associated with a charged lepton and two quarks, which can induce both production and decay of the heavy neutral leptons (HNLs) simultaneously. For HNLs of O(10) GeV mass, such operators with a NP scale above ∼ 1 TeV can easily make the HNLs become long-lived, leading to displaced vertices at the LHC. We have therefore proposed a displaced-vertex search strategy based on a prompt-lepton trigger and selection of highquality displaced tracks. By performing Monte-Carlo simulations with MadGraph5 and Pythia8, we have estimated the sensitivity reaches for ATLAS in the high-luminosity LHC era with 3 ab −1 integrated luminosity, to four single-N R EFT operators: O duN e , O LN Qd , O LdQN , and O QuN L .
Multiple combinations of quark and lepton flavors can be studied. Here, we have considered mainly two combinations: (u, d) and (c, s). The first (second)-generation-quark only flavor combination is then projected to have the best (worst) sensitivities, because of their portion in the proton content. For both quark combinations, we also studied all possible lepton generations, i.e. electron, muon and tau. In addition, for simplicity, we did not take into account the effect of the active-sterile neutrino mixing, which is supposed to be negligible if the type-I seesaw relation is assumed. For the (u, d) and (c, s) combinations, we find in general for the considered single-N R operators, ATLAS can probe Λ up to 20 TeV and above for m N 20 GeV, if we switch on one operator at a time.
In addition to the EFT scenarios, we also revisited the minimal scenario of the HNL mixing with the SM neutrinos. In this scenario, the type-I seesaw relation is not assumed and we have two independent parameters: mass of the HNL and its mixing parameter with one type of the active neutrinos: a simple 3+1 scenario. These results are an update of those given in ref. [52]. Besides some minor changes, the most important difference is that we have now taken into account both charged and neutral currents in our computation, leading to more realistic projection results, especially for the case of mixing with the τ neutrino.
In summary, we conclude that a displaced-vertex search at ATLAS for HNLs can probe new physics scales up to about 20 TeV and, in some cases above, for HNL mass between about 5 GeV and 50 GeV, depending on the quark and lepton flavors associated with the single-N R operator under consideration.