Heavy neutral fermions at the high-luminosity LHC

Long-lived light particles (LLLPs) appear in many extensions of the standard model. LLLPs are usually motivated by the observed small neutrino masses, by dark matter or both. Typical examples for fermionic LLLPs (a.k.a. heavy neutral fermions, HNFs) are sterile neutrinos or the lightest neutralino in R-parity violating supersymmetry. The high luminosity LHC is expected to deliver up to 3/ab of data. Searches for LLLPs in dedicated experiments at the LHC could then probe the parameter space of LLLP models with unprecedented sensitivity. Here, we compare the prospects of several recent experimental proposals, FASER, CODEX-b and MATHUSLA, to search for HNFs and discuss their relative merits.


Introduction
Light long-lived particles (LLLPs) appear in many extensions of the standard model (SM). LLLPs can be scalars, fermions or vectors. Fermionic LLLPs are also often called heavy neutral fermions (HNF) in the literature. LLLPs are usually motivated either by dark matter, by small neutrino masses or by both. A discussion of different theoretical models for LLLPs can be found, for example in [1,2].
In the past few years several experimental proposals with improved sensitivities to LLLPs have been discussed. For example, there are the planned fixed target experiment SHiP [2], the near detector of the future DUNE experiment [3], or also NA62 [4]. However, note that the primary goal of NA62 is to measure precisely Br(K + → π + νν), while the main task of the near detector of DUNE is just monitoring the neutrino flux for the far detector [3].
It is expected that the LHC will deliver up to L = 3000/fb of luminosity over the next (15)(16)(17)(18)(19)(20) years [5]. Perhaps unsurprisingly, a number of new proposals to search for LLLPs have appeared, all based on the idea to exploit LHC's large luminosity: MATHUSLA [6], CODEX-b [7] and FASER [8]. The physics potential of these three experiments has so far not been fully discussed in the literature and it is the aim of the current paper to estimate, and compare to each other and previous experiments, the sensitivity of these proposals for fermionic LLLPs.
MATHUSLA [6] is a proposed very massive detector, possibly to be located above ground on top of the ATLAS experiment. The sizeable distance of MATHUSLA from the interaction point (IP) of the LHC beams will allow to test for rather long life-times, up to cτ < ∼ 10 (7−8) m, but requires a huge detection volume. CODEX-b [7] is a proposal that takes advantage of a relatively large shielded empty space near the location of the LHCb experiment. Being closer to the IP, CODEX-b proposed size is much smaller than that of MATHUSLA. Finally, the authors of FASER [8] propose to construct a very modestly sized detector, situated in the very forward direction, close to either the ATLAS or the CMS IP. Note that FASER is discussed in three different variants in [8,9]. We will summarize the experimental parameters of the different FASER setups and those of CODEX-b and MATHUSLA in section 3.
In the original paper on the MATHUSLA detector [6], exotic Higgs decays to LLLPs were used to study the sensitivity of the experiment. Similarly, the authors of CODEX-b [7] used a Higgs portal model for LLLPs to investigate the reach of their proposal. They also considered [7] a light neutral scalar, that mixes with the Higgs, produced in B-mesons decays as a LLLP candidate. The original FASER [8] publication studied dark photons produced in meson decays to estimate the sensitivity of the different FASER setups. Here, we study the reach of these three experimental proposals for the case that the LLLP is a heavy neutral fermion (HNF). We concentrate on two particular example models of HNFs: (i) Sterile neutrinos and (ii) the lightest neutralino in R-parity violating supersymmetry. As discussed in section 2 both models are motivated by being possible explanations for the observed small neutrino masses (and mixings). In our simulation we consider all possible LLLP production channels: D-mesons, B-mesons, W and Z bosons, as well as Higgs boson. We compare the different experiments systematically for the different channels and then discuss sensitivities for our example models.
We note that, very recently in [10] the sensitivity of FASER to sterile neutrinos was estimated. We will comment on this work in more detail later, but note here only briefly that our estimates roughly, but not completely, agree with those given in [10]. However, [10] studies only the case of sterile neutrinos, while we also discuss R-parity violating neutralinos, and concentrates exclusively on FASER, while we take into account the different experimental proposals discussed above and also compare to the beam-dump experiments [2][3][4].
Before closing this introduction, we mention that there exist of course already many searches for sterile neutrinos (and other HNFs). For a review on constraints for sterile neutrino see, for example [11]. Also the main LHC experiments, ATLAS and CMS, have searched for HNFs. ATLAS [12] published results of a search based on the final state lljj, giving only weak upper limits on the mixing of the sterile neutrinos V 2 αN (10 −2 − 10 −1 ) (α = e, µ) for m N ∼ (100 − 500) GeV. CMS searched for sterile neutrinos in trilepton final states and very recently published limits as low as V 2 αN 10 −5 in the mass range (10−100) GeV [13]. With these results [13], CMS now gives limits competitive with those derived by the DELPHI experiment at LEP [14].
Note that for small mixing angles V 2 αN below, say V 2 αN ∼ 10 −7 , for m N ∼ O(10) GeV, the decay lengths of sterile neutrinos become large enough to be detected exerpimentally and ATLAS/CMS could search for sterile neutrinos using the "displaced vertex" signal [15]. However, current displaced vertex search strategies, as used by CMS [16] for example, are not very well suited for light, say m N < ∼ 100 GeV, sterile neutrinos [17].
We also mention in passing that our other HNF candidate, the neutralino, has been studied as an LLLP candidate before. R-parity violating SUSY and neutralinos as LLLPs are mentioned in the SHiP proposal [2] and the SHiP sensitivity for neutralinos has been studied in more details in [18,19].
The rest of this paper is organized as follows. In section 2 we briefly summarize the basics of sterile neutrinos and neutralinos with R-parity violation. In section 3 we give a basic description of the different experiments and outline our simulations. In section 4 we discuss our numerical results, before closing with a short summary.

Heavy neutral fermions: Example models
In this section we give a short summary of the two models for heavy neutral fermions that we study numerically in this paper. We first discuss sterile neutrinos and then give some basic definitions and features for the neutralino in R-parity violating supersymmetry. Since both models have been discussed in the literature many times, we will be very brief.

Sterile neutrinos
The standard model predicts neutrinos to be massless, in contrast to the results of neutrino oscillation experiments. 1 The simplest extension of the SM, which can explain the experimental data, adds n fermionic singlets. Oscillation data requires n ≥ 2. The Lagrangian of this model contains two new terms Here, we have suppressed generation indices. In general M N is a complex symmetric (n, n) matrix, while Y ν is a (3, n) matrix. In the simple model considered here, without new interactions for the ν R , one can perform a basis change and choose the entries of M N to be diagonal, real and positive. The masses of the active neutrinos are small, if (Y ν v)·M −1 N 1, this is the essence of the seesaw mechanism. Diagonalization of the mass matrix leads then to three light, active neutrinos and n nearly sterile mass eigenstates, which we denote by ν S in the following.
The heavy sterile neutrino Charged (CC) and Neutral Current (NC) interactions are where i = 1, 2, 3 and j = 1, .., n and α denotes the charged lepton generation. The lefthanded sector neutrino mixing matrix V L is measured in neutrino oscillations. V αN j describes the mixing between ordinary and sterile neutrinos. Within the simple seesaw model, described by eq. (2.1), one expects that V αN j is roughly of the order of V αN j ∝ m ν /M N , i.e. |V αN j | 2 5 × 10 −11 ( mν 0.05eV )( 1GeV M N ). However, in extensions of this simple framework, for example the inverse seesaw [21], much larger values for the mixing can occurr, despite the smallness of the observed neutrino masses. For this reason, for the sensitivity estimates of the different experiments we will take |V αN j | 2 as a free parameter in our calculations. Note that the mixing between sterile and active neutrinos controls both, production and decay of the sterile states.
Oscillation data shows two large mixing angles in the active neutrino sector [20]. Thus, one expects that the heavy sterile neutrinos couple typically to more than one generation of charged leptons too, see eq. (2.2). It is easy to fit all oscillation data with the seesaw mechanism, described by eq. (2.1). However, the Yukawa matrices can be fixed by such a fit only up to an orthogonal rotation matrix containing three complex parameters [22], leaving |V αN j | essentially as free parameters. 2 In our sensitivity estimates we will simply assume that only one sterile neutrino exists in the mass range to which the experiments are sensitive. We will also not distinguish between e and µ flavours, assuming simply that only one of the corresponding |V αN | is non-zero. Since we are only interested in estimating sensitivity ranges, not in a full reconstruction of the seesaw parameters, this should be a reasonable approximation.

R-parity violating neutralino
Supersymmetric models with R-parity violation (RPV) have been discussed in detail in many publications in the literature, for reviews see for example [24,25]. R-parity violating terms do either break baryon (B) or lepton (L) number. The lepton number violating (LNV) part of the superpotental can be written as: while the baryon number violating (BNV) terms are: It is experimentally excluded that both terms are present at the same time, since otherwise the proton decays with an unacceptable rate unless the product of the two couplings is tiny [25], typically λ × λ < ∼ 10 −24 for TeV SUSY masses. The lepton number violating terms generate neutrino masses and it is well-known that even the simplest bilinear RPV model can fit oscillation data [26]. We will therefore put all BNV terms to zero in the following. Once R-parity is violated, the lightest SUSY particle is no longer stable and therefore there are no constraints on its nature from cosmology. Thus, even charged or coloured SUSY particles could be the LSP (lightest supersymmetric particle). Here, we are exclusively interested in the case where the lightest neutralino is the LSP.
In the so-called CMSSM the gaugino mass terms M 1 and M 2 follow the approximate relation M 1 (1/2)M 2 . This leads to a lower limit on m χ 0 1 of roughly m χ 0 1 > ∼ 46 GeV [27]. 2 For an extension of this Casas-Ibarra parametrization for the inverse seesaw case, see [23].
However, in more general SUSY models, M 1 and M 2 are just free parameters and it is easy to show that for [28] the lightest neutralino is massless at tree-level. In our numerical studies we will simply take the mass of the lightest neutralino, m χ 0 1 as a free parameter, without resorting to any underlying SUSY breaking model. Note, however, that this lightest neutralino necessarily has to be mostly bino, due to the lower mass limits on charginos that can be derived from LEP data [27].
Decays of the lightest neutralino can be induced via either bilinear or trilinear RPV terms. In the former case, the neutralinos and the neutrinos mix at tree-level, thus the neutralino can decay via diagrams involving W and Z bosons. A rough guess for the width of the neutralino in BRPV (bilinear RPV) can be derived from the results of [29]. An order of magnitude estimate can be given as: Contributions to the decay of the neutralino from trilinear RPV terms, on the other hand, involve sfermion exchange diagrams. The order of magnitude of the decay width of the neutralino can be estimated in this case as: Here, λ stands symbolically for any of the trilinear couplings in eq. (2.3) and mf is the corresponding scalar fermion mass (either squark or slepton). Given that there are currently only lower limits on all sfermion masses, it is possible to take either trilinear couplings and the sfermion masses or simply the total width of the χ 0 1 as a free parameter in trilinear RPV.
At the LHC, neutralinos can either be pair produced, via diagrams involving either a Z-boson or a Higgs, or singly produced via R-parity violating couplings. Since one expects R-parity violating couplings to be small parameters [24][25][26], we will focus on pair production. Z bosons are produced much more abundantly in the LHC than the Higgs. We will therefore concentrate the following discussion on Z-bosons. Z bosons can decay to pairs of neutralinos, if m χ 0 1 < ∼ m Z /2. The decay width Γ(Z → χ 0 1 χ 0 1 ) has been calculated in [30]. Important for us is the coupling between Z-boson and two neutralinos: Here, N ij is the matrix that diagonalizes the neutralino mass matrix and c 2β /s 2β are cosine and sine of β, with tan β being the usual ratio of the vacuum expectation values of the two Higgs doublets. Thus, the relevant coupling for the decay to neutralinos is proportional to the Higgsino content in χ 0 1 and is not suppressed by small RPV parameters. Different from the case of the sterile neutrino, therefore for neutralinos production cross section and decay length are not related. This important distinction between our two LLLP candidates will be discussed in more detail in section 4.3.
Importantly, there is an upper bound on Γ(Z → χ 0 1 χ 0 1 ) from the LEP measurement of the invisible width of the Z-boson. The PDG [27] gives for Γ(Z →inv) a value in agreement with the standard model calculation with three generations of light neutrinos and the error bar on the measurement corresponds to an upper limit on the branching ratio into additional invisibly final states of roughly Br(Z → χ 0 i χ 0 j ) < ∼ 0.1 % at 90 % c.l. The Higgsino content in the lightest neutralino depends mostly on the parameter combination M 1 /µ. LEP gives a lower limit of roughly µ > ∼ 100 GeV [27]. A recent LHC search in ATLAS for electroweak SUSY production [31] excludes now values of µ, depending on other SUSY parameters, up to roughly µ 130 GeV. For a cross check, we used the model MSSMTriRpV from the repository of SARAH-4.12.2 [32]. We perform numerical calculations with SPheno-4.0.3 [33,34]. For the choice M 2 = 500 GeV, µ = 130 GeV, tan β = 10, (2.9) and a lightest neutralino mass m χ 0 1 m Z /2, we find Br(Z → χ 0 i χ 0 j ) 0.06 %. Thus, given current constraints on SUSY parameters, the Higgsino content in the lightest neutralino can still be large enough to (nearly) saturate the experimental bound on Br(Z → χ 0 i χ 0 j ). 3 In our numerical calculations, see section 4.3, we will not do a scan over the soft SUSY breaking parameters. Instead we will treat both, the mass of the lightest neutralino and Br(Z → χ 0 1 χ 0 1 ) as a free parameter in our numerical study. Note, however, that a future lower limit on µ larger than the numbers quoted above will consequently result in smaller values for the maximally achievable Br(Z → χ 0 1 χ 0 1 ). We would also like to mention that the lower limits on charged SUSY particles require that the lightest neutralino must be mostly bino, for the low mass we consider in section 4.3.

Experimental setups & simulation
Here we first give a brief description of the different experiments considered. For more details we refer to the original publications for MATHUSLA [6], CODEX-b [7] and FASER [8]. We then describe our numerical simulation of these proposals. Note that, given that the experimental setups might still undergo some changes and refinements we have not aimed at very high accuracy in our simulations of the experiments, although we have checkedwherever possible -that our calculations agree with results obtained in previous works, see also section 4.
CODEX-b ("Compact detector for Exotics at LHCb") [7] was proposed recently as a detector for LLLPs. The experiment consists of a cubic box with approximate dimensions of (10 × 10 × 10) m, installed in a free space near the LHCb experiment. As discussed in [7], with some modest amount of additional shielding CODEX-b is expected to operate in the low background environment, necessary to search for very rare events.
The FASER ("ForwArd Search ExpeRiment") proposal uses a cylindrical detector a few hundred meters downstream of the ATLAS or CMS IP. The FASER papers [8][9][10] discuss several different options for the position and size of the detector. We use the following three options given in [9]: FASER r , the small FASER with radius 0.2 m, FASER R (large Radius) with radius of 1 m and FASER n ("near") with small radius 0.04 m at a shorter distance. Note that especially for FASER n , backgrounds from the LHC beam might be a concern, since there is rather little shielding in the FASER near position.
MATHUSLA ("MAssive Timing Hodoscope for Ultra Stable neutraL pArticles") [6] proposes to put a surface detector above the ATLAS interaction point. With a distance between 140 − 320 m from the IP, the MATHUSLA detector has to be quite massive, compared to the other two proposals. Its dimensions are given as [6] 200m × 200m × 20m. Being above ground, cosmic rays are a serious background concern. The authors of [6] discuss in some detail, how this background can be kept under control, surrounding the decay volume with scintillator detectors and Resistive Plate Chambers (RPCs). The excellent timing resolution of these anti-coincidence detectors will then allow to cut cosmic ray backgrounds to near negligible levels.
In Table. 1, we summarize the relevant parameters of the different experiments. L is the value of the luminosity assumed in the original proposals. Note that the authors of [7] give an expected luminosity of 300/fb for LHCb for the end of Run-5 of the LHC. 4 This is a factor of 10 lower than used in the other proposals, which make use of the high luminosity environment in ATLAS or CMS. We also give f D , the maximum fraction of events decaying within the detector length, estimated for the optimal decay length, given the corresponding L min and L max . geometric is the geometric factor, calculated relative to the full solid angle. Even FASER R corresponds to an instumented volume of only roughly 8 m 3 , while CODEXb corresponds to 10 3 m 3 and the massive MATHUSLA covers a volume of 8 · 10 5 m 3 . Note, however, that FASER sits at very large η. Thus, geometric underestimates the sensitivity of FASER for events coming from D-and B-meson decays, see the discussion in the next section. We note that in [10], the authors discuss a slightly modified setup for FASER with respect to those quoted in table (1). In particular they use a distance of L max =480 m to the IP. We will comment on the resulting differences with respect to our calculation in section 4.2.
We now briefly describe our simulation. We use SARAH-4.12.2 [32] to generate UFO [39] models for sterile neutrinos in a Seesaw Type-I model and the RPV-MSSM. We generate spectrum files with SPheno-4.0.3 [33,34]. We created the Seesaw Type-I model using the SARAH environment, while the RPV-MSSM model already exists in the model repository.
For the W, Z and Higgs boson channels, we import Seesaw Type-I or RPV-MSSM UFO models into MADGRAPH5 aMC@NLO v2.6.0 [40] where these bosons are generated, put on-shell and decay. We then read the LHE [41] event records with MadAnalysis5 v1.5 [42][43][44], which we use to apply the relevant geometric cuts on the pseudo-rapidity. With these simulated events we calculated the number of LLLPs produced in the detector window (n f ) and the mean values of their βγ. ( β · γ is needed for the correct simulation of the decay length.) Dividing n f by the total number of simulated events provides f window , the fraction of events relative to the total cross section. We prefer to use this procedure, since we expect that the values of the total cross section simulated in MadGraph5 should have a larger uncertainty than the relative fractions, as obtained in our procedure.
For the meson channels, we use the HardQCD:hardccbar( HardQCD:hardbbbar) matrix element calculator with Pythia 8.205 [45,46] for generating events gg, qq → cc(bb), showering and hadronization. For simulating sterile neutrinos, we need to know the fraction (f window ) of events of the total cross section moving towards the different detectors and the mean β · γ of these windowed sterile neutrinos. These can all be calculated and given by Pythia8. In particular, Pythia8 allows to add a fourth neutrino to the standard model module, but does not automatically recalculate branching ratios, when one varies the mass of this fourth neutrino. For our sensitivity estimates we thus include the phase space suppression of these branching ratios by hand, which is important for non-zero sterile neutrino masses. We allow the B-and D-mesons to three-body decay to the fourth neutrino plus an electron and a lighter meson. Moreover, for a sterile neutrino mass close to that of B-or D-mesons, the leptonic two-body decay of the charged pseudoscalar mesons (in our case, D ± , D ± s , B ± , and B ± c ) can lead to relatively large contritbution. We therefore also include them in our calculation.
We then calculate the total number of events, using mostly experimental data. For the D-meson and B-meson production we use the results from LHCb [47-50]. LHCb gives cross sections within certain ranges of pseudo-rapidity η and p T . We therefore use FONLL [51][52][53][54] to extrapolate these cross sections to the full range of η and p T . We checked that FONLL estimates roughly an uncertainty order 15 % for the total cross section in case of D-mesons (and a similar error for B-mesons).
Ref. [55] gives cross section of the W ± → l ± ν and Z → l + l − production (where l ± = e ± , µ ± ) in proton-proton collisions at √ s = 13 TeV. Making use of the branching ratios listed in the PDG [27], we obtain the total cross sections for W and Z bosons production (and also for the D-and B-mesons). The Higgs boson is produced at the LHC dominantly through gluon fusion. The cross section can not be measured independently from the Higgs decay branching ratios, thus the most reliable value for the total cross section is probably still the calculated value of 43.92 pb, as given in [56]. With the total cross sections fixed this way, we need only f window and β · γ from our simulation to calculate the total number of events in the different detectors.
We calculate the decay width of the sterile neutrinos according to formulas given in [11]. This should give a more accurate treatment than the decay width calculated in SPheno, which does not take into account hadronic form factors for decays to mesons. We have checked that these form factors are numerically important for sterile neutrino masses below 5 GeV. For the decay width of the lightest neutralinos we take directly the output from the SPheno spectrum files. 5 We calculate the number of LLLPs decaying inside the detector using the following formula, (3.1) Here, N decay is number of the LLLPs decaying in the detector, N total is the total number of LLLPs produced at the IP, while L decay = βγcτ is their decay length. f window is the fraction of events of the total emitted into the detector window. L min (L max ) is the minimum (maximum) distance from the IP to the detector. For each detector, N total and f window are functions of the mass of the LLLPs, and as described above, we simulate these values using either MadAnalysis5 or Pythia8.
The simple description given in eq. (3.1) corresponds to approximating the detectors as cones defined by their coverage in η and φ. Since the orientation of MATHUSLA is not well described by such a simple cone with a tip at the IP, we refine this part of the calculation for MATHUSLA. We assume that f window and β ·γ are constant over the range of η covered and then divide MATHUSLA into 10 smaller boxes of equal size. We then sum up the decay number of events in each box to obtain the total number. We have tested this rather simple approximation against a calculation done with a Mathematica notebook, distributed by the authors of [6], which integrates numerically over the detector valume. We find good agreement between our simple-minded approach and the more accurate one given in [6].

Sensitivities for different production modes
Sterile neutrinos are produced via their mixing with the active neutrinos. Thus, for sterile neutrinos production and decay width are both governed by the same parameter, ie. |V αN j | 2 . The situation is different for other LLLP candidates. For example, as discussed in section 2.2, the neutralino in RPV SUSY can be singly produced via its mixing with the neutrino or produced in pairs, via Z-or Higgs diagrams. Since one expects that RPV parameters are small (since neutrino masses are small) also neutralino-neutrino mixing is expected to be small. Thus, the main production mode for a light neutralino should be Z-boson (and/or Higgs) boson decays.
For this reason, we will first give sensitivity estimates for the different experiments separated into the different possible production modes. We consider D-mesons, B-mesons, the gauge bosons W and Z, as well as the Higgs. Also, in this subsection we will present our estimates in the particular parameter plane branching ratio versus cτ . This has the advantage that the results shown can be easily used to estimate also sensitivities for other LLLP candidates.
In fig. (1) we show the estimated reach of CODEX-b, FASER and MATHUSLA in the plane branching ratio versus cτ for D-mesons (left) and B-mesons (right). Here, and in all other plots in this paper we use the contours of 4 signal events for estimating the sensitivity limit. 6 As the plot shows, FASER R and MATHUSLA are sensitive to different regions in cτ , despite being at similar distances from the IP. This can be understood, because mesons flying in the foward direction receive typically much larger boosts than mesons produced more centrally. The same effect makes FASER less sensitive at large cτ . One notes that at cτ 1 m FASER R can reach nearly as small branching ratios as MATHUSLA at cτ 10 2 m in the case of D-mesons, despite being a much smaller experiment. The reason is simply that charm (and to a lesser degree bottom) quark production at the LHC is very strongly peaked in the forward direction. For the same reason the different setups of the FASER experiment (for parameters see table (1)) do particularly well for D-mesons: while FASER r contains only 4 % of the decay volume of FASER R , its sensitivity is only a factor of 6 smaller for events coming from D-meson decays. FASER n is sensitive to smaller cτ , but is always less sensitive than the other FASER variants. CODEX-b is less sensitive than either FASER R or MATHUSLA in case of meson decays. However, this is partly due to the lower luminosity (300/fb versus 3000/fb) assumed for CODEX-b. Finally note that MATHUSLA does significantly better than FASER R for B-mesons. In fig. (2) we show the results for production from W-(left) and Z-bosons (right). Note the change of scale in the axis: This simply reflects that there are 5 orders of magnitude more D-mesons than gauge bosons produced at the LHC. Comparing the results shown in fig. (2) with those of fig. (1) one sees that CODEX-b now does better than even FASER R for both, W and Z-bosons. MATHUSLA is again the most sensitive setup, apart from some region of parameter space at small cτ . Note also that FASER r and FASER n give much weaker limits than FASER R in this case. The explanation is again simply that for gauge bosons there is much less boost into the forward direction than for D-mesons.
We have also calculated the corresponding expectation for the different experimental  setups for Higgs production, see fig. (3). The figure shows no contours for FASER r and FASER n since these variants have no sensitivity in the plane plotted here (even with 3000/fb of luminosity). Also FASER R can not compete with CODEX-b or MATHUSLA in case of Higgs production. CODEX-b is around a factor ∼ 40 less sensitive than MATHUSLA. However, recall that (i) CODEX-b is a much smaller setup and (ii) this estimate uses only 300/fb for CODEX-b. In fig. (3) we show the sensitivity ranges for two choices of m νs . To the left we use m νs = 1 GeV, while the plot on the right uses m νs = 10 GeV. This change of mass leads to a change in the value of β · γ and thus to a corresponding shift in the region in cτ , where the experiments are most sensitive. Since this shift is similar for all the three experiments shown, however, this does not affect the conclusions. We note in passing that the curve for MATHUSLA on the right of fig. (3) is very similar to the corresponding one shown in fig. 3 of [6] for a scalar LLLP with the same mass.

Sterile neutrinos
We now turn to a discussion of the sterile neutrino results. For the calculation of the sensitivities of the different experimental proposals we take into account the different production processes discussed in the previous subsection: D-meson and B-meson decays, Wand Z-bosons, as well as Higgses. For the sterile neutrino case, production of steriles from Higgs decays gives only a negligible contribution to the sensitivity, as expected.
In our estimates we will not consider sterile neutrinos decaying to τ 's. We will also consider only mixing of the sterile neutrinos with either e or µ for simplicity. For this simplified case, assuming the experimental detection efficiencies for e and µ to be similar, plots for |V eN | 2 and |V µN | 2 are very similar and we simply will show plots using |V αN | 2 .
In fig. (4) we show a comparison of the different experiments in the plane |V αN | 2 versus mass of the sterile neutrino, m ν S [GeV]. To the left we compare the different variants of FASER to each other, while the plot on the right compares FASER R , CODEX-b and MATHUSLA. We also show the projected sensitivity of SHiP [2,57], the expectations for the near detector of the future DUNE experiment [3], denoted as LBNE in the plot, and a recent estimate for the final sensitivity of NA62 [58]. For a description of the NA62 experiment see [4]; for the current status of limits for HNFs from NA62 see [59]. The grey area in the background is excluded from past searches [60]. The experimental references used in drawing the excluded areas are PS191 [61], JINR [62], CHARM [63] and DELPHI [14].
The plot on the left of fig. (4) shows that, as expected, FASER R always does better than the other FASER variants. It also shows that below m ν S 2 GeV FASER R is competitive with NA62 [4], while for m ν S > ∼ 2 GeV, FASER has better sensitivity than NA62. For the whole mass range, SHiP is more sensitive than FASER R .
The plot on the right of fig. (4) shows that CODEX-b and FASER R have quite similar sensitivities below m ν S 3.2 GeV to sterile neutrino parameters, while MATHUSLA does better than both in most parts of the parameter space. While the plot shows that SHiP has the best expected sensitivity in the mass range 0.5 < ∼ m ν S < ∼ 2 GeV, for larger masses FASER R and CODEX-b are only worse than SHiP by approximately one order of magnitude. One notes also that MATHUSLA is only slightly less sensitive than SHiP for m ν S < ∼ 2 GeV, and even better than SHiP for 2 < ∼ m ν S < ∼ 4 GeV. Below m ν S < ∼ 0.5 GeV the most sensitive experiment will be LBNE [3]. Note, however, that our calculation does not include sterile neutrinos from Kaon decays, which provide most of the sensitivity of LBNE.
We also briefly compare the results of [10] with our calculations. As mentioned in the introduction, [10] estimates the sensitivity of FASER for sterile neutrinos. While our results are overall similar to those shown in [10], there are also some minor differences. We believe this is mainly due to the following: (i) [10] gives a decay length around a factor of 2 larger than our calculation in this mass range; (ii) the distance of FASER in [10] is set to 480 m (we use 400 m from the original proposal).

The neutralino case
We now turn to a discussion of light neutralinos in R-parity violating SUSY as LLLP candidates. In SUSY models with BRPV, neutralinos and neutrinos mix. As discussed in section 2.2, this mixing is related to the neutrino mass generated in these models. Constraints on purely BRPV models for singly produced neutralinos can thus be derived from a re-interpretation of the sterile neutrino constraints discussed above.
The main phenomenological difference between sterile neutrinos and neutralinos in RPV then actually comes from pair production of neutralinos. As discussed in section 2.2, pair production of the lightest neutralino is not suppressed by small RPV parameters. We will consider neutralinos that are pair-produced from Z-boson decays. At the LHC with 3000/fb of luminosity there will be a total of more than 10 11 Z-bosons produced. The current upper limit of Br(Z → χ 0 i χ 0 j ) < ∼ 0.1 % at 90 % c.l. implies that at the maximum allowed still up to 10 8 neutralinos could be produced from Z-decays.  The plot on the left puts the branching ratio at the experimental upper limit, while the plot on the right uses Br(Z → χ 0 1 χ 0 1 ) = 10 −5 . 7 As the figure shows, for Br(Z → χ 0 1 χ 0 1 ) = 10 −5 the sensitivity of FASER R gets significantly diminished. Roughtly for Br(Z → χ 0 1 χ 0 1 ) = 5 · 10 −6 there will be less than 4 events in FASER R for any combination of cτ and m χ 0 1 from pair produced neutralinos. For FASER r and FASER n the corresponding limits are 3 · 10 −4 and 10 −4 , respectively. Thus, the smaller variants of FASER could have sensitivity only, if Br(Z → χ 0 1 χ 0 1 ) is very close to the current upper bound. For this reason, we do not show contours for these two variants in fig. (5). For CODEX-b and MATHUSLA the numbers are much more optimistic; we estimate these two experiments can probe values down to Br(Z → χ 0 1 χ 0 1 ) = 6 · 10 −7 and 1.4 · 10 −8 respectively.
As fig. (5) also shows the parameter space testable in FASER R , CODEX-b and MATH-USLA covers the range of decay lengths expected for such light neutralinos in SUSY models with bilinear RPV, see also the discussion below eq. (2.6) in section 2.2. It is important to note that for pair produced neutralinos, a much larger mass range can be covered than in the sterile neutrino case. As the figure shows, neutralinos up to m Z /2 can be tested. Again, this is due to the fact that production and decay of the neutralino are not related in pair production. Note also that, in particular, MATHUSLA can probe large values of cτ , which are interesting if (i) the neutralino is very light and/or (ii) the SUSY parameters M 2 and/or µ are large. Recall, however, that large values of µ automatically imply a small higgsino content in χ 0 1 , leading to small values for Br(Z → χ 0 1 χ 0 1 ). We now turn to a discussion of the case of trilinear RPV. In fig. (6) we compare the different experiments in the plane λ ijk versus lightest neutralino mass for different values of other parameters. Here, for λ ijk the generation indices i, j, k could take in principle any value 1, 2, 3, depending on the final state discovered. However, in practice the results shown are valid only for the first two generations of quarks and leptons, since we do not take into account phase space suppression due to non-zero final state quark and lepton masses in our calculation. The top row shows the sensitivities for two different values of Br(Z → χ 0 1 χ 0 1 ). Here, we have fixed all scalar masses (both sleptons and squarks) to m SU SY 2 TeV. The two plots in the bottom row show how the explorable regions change for larger values of the sfermion masses. In these plots we assume that the decay length is dominated by trilinear RPV and neglect any possible contribution from BRPV.
Since the smallness of λ ijk controls only the decay length, but not the production cross section, in principle very small values of λ ijk are accessible in these searches. Note that the values of λ ijk shown can reach values several order of magnitudes smaller than even the best of the current current upper limits on λ ijk [24,25]. However, recall that pair production depends on the unknown value of Br(Z → χ 0 1 χ 0 1 ) , as discussed above. Again, FASER R is expected to be less sensitive than CODEX-b, with the best limits expected from MATHUSLA.
In summary, if the lightest neutralino has a mass in the range (few) GeV to m Z /2, LLLP searches at FASER, CODEX-b and MATHUSLA can probe part of RPV SUSY parameter space not accesible in any other experiment. In particular, for bilinear RPV MATHUSLA can cover large part of the predicted range of decay lengths for such light neutralinos.

Conclusions
We have discussed the sensitivities for three recent experimental proposals, MATHUSLA [6], CODEX-b [7] and FASER [8] for the case of fermionic light long-lived particles. We considered two concrete example models for our study: (a) light sterile neutrinos and (b) the lightest neutralino in R-parity violating supersymmetry. Both candidates are motivated by theoretical models that can explain the observed small neutrino masses.
For sterile neutrinos, FASER R and CODEX-b show similar sensitivities. Here, FASER R compensates its smaller detection volume by taking advantage of the fact that D-mesons (and to some degree B-mesons) are produced mostly in the forward direction. MATHUSLA is more sensitive than either FASER R or CODEX-b and, in fact, is competitive with the fixed target experiment SHiP.
For the case of neutralinos in RPV SUSY we have found that FASER R , CODEX-b and MATHUSLA can cover interesting parts of the parameter space of these models, if the lightest neutralino has a mass in the range of a few GeV up to m Z /2. In particular, for bilinear RPV models these experiments cover large parts of the range of cτ predicted theoretically from the observed neutrino masses. For trilinear RPV, on the other hand, we have shown that if such a light neutralino exists, RPV couplings can be probed which are orders of magnitude smaller than all existing constraints.
Quite generally, MATHUSLA [6] shows better sensitivity to fermionic LLLPs than either CODEX-b [7] or FASER [8]. However, this advantage clearly comes at a price: MATHUSLA has by far the largest instrumented volume of all three experiments. Considering that the variants of FASER, discussed so far in the literature [8,10] are actually quite small, compared to the other experiments, we think it would be very interesting to study, if space for a larger version of FASER in the very forward direction could be found.