Heavy Neutral Leptons at FASER

We study the prospects for discovering heavy neutral leptons at ForwArd Search ExpeRiment (FASER), the newly proposed detector at the LHC. Previous studies showed that a relatively small detector with ~10 m length and ~1 m cross sectional area can probe large unconstrained parts of parameter space for dark photons and dark Higgs bosons. In this work we show that FASER will also be sensitive to heavy neutral leptons that have mixing angles with the active neutrinos that are up to an order of magnitude lower than current bounds. In particular, this is true for heavy neutral leptons produced dominantly in $B$-meson decays, in which case FASER's discovery potential is comparable to the proposed SHiP detector. We also illustrate how the search for heavy neutral leptons at FASER will be complementary to ongoing searches in high-$p_T$ experiments at the LHC and can shed light on the nature of dark matter and the process of baryogenesis in the early Universe.


I. INTRODUCTION
In the past years, the Large Hadron Collider (LHC) has collected an impressive amount of data and placed constraints on multitude of models for new physics. Most of these searches are targeting high-p T signatures corresponding to new strongly interacting heavy particles. However, up to now there has been no discovery of any elementary beyond the standard model (BSM) particle, which has motivated the community to consider a broader range of different physics signatures. In particular, if new particles are light and weakly coupled, they could travel a macroscopic distance before decaying. Searches for such long-lived particles (LLPs) are experimentally clean and theoretically well-motivated (see, e.g., [1]). While current and planned searches for MeV to GeV range LLPs typically employ beam-dump experiments or meson factories, LLPs would also be produced abundantly in high-energy pp-collisions at the LHC. However, they are typically produced with low-p T and therefore will move along the beam pipe and escape detection in the ATLAS [2] and CMS [3] experiments.
In recent papers [4,5], we proposed a new experiment searching for LLPs, FASER, the ForwArd Search ExpeRiment at the LHC. FASER would be placed in the forward direction from the ATLAS or CMS interaction point (IP) and operate concurrently with the ongoing high-p T searches. In particular, we considered a small size detector (20 cm to 1 m radius, 10 m long cylinder) placed in a representative location 400 m away from the IP along the beam axis, possibly in a side tunnel, after the main LHC tunnel enters the arc section and starts to curve. Importantly, the existing LHC infrastructure, including neutral absorbers and magnets deflecting charged particles, would provide natural shielding from various standard model (SM) backgrounds. An additional layer of shielding is provided by the rock and concrete that separate the main tunnel from the side tunnel in which FASER would be placed.
A more detailed study of the LHC infrastructure around the ATLAS IP revealed an attractive location for FASER in the currently unused side tunnel TI18 in the distance about 480 m away from the IP. It is a former service tunnel used by LEP as a connection between the SPS and the main tunnel. In the following, we will perform sensitivity studies for this location.
In this study, we analyze FASER's prospects for discovering heavy neutral leptons (HNLs) as a well-known example of new fermionic particles that can be found in forward searches [6][7][8][9]. This analysis is complementary to our previous studies. In particular, unlike dark photons, HNLs with mass above ∼ 1 GeV can be abundantly produced in heavy meson decays thanks to their mixing with the SM (active) neutrinos. However, in contrast to dark Higgs bosons, HNLs couplings to mesons are not dictated by the quark Yukawa couplings, but rather by the set of unknown Yukawa couplings between HNLs and the active neutrinos, as well as by the CKM mixing parameters. As a result, HNLs can also be efficiently produced in D meson decays, while less abundantly produced B mesons play a dominant role only when other channels are kinematically forbidden.
Alternatively, one can also consider a detector in some near location inside the straight intersection of the main LHC tunnel between the two beam pipes after they split. The location that is as close as possible to the IP is right behind the TAN neutral particle absorber which we chose as a representative case in [4,5]. The expected signal of new physics in this case can be enhanced with respect to the far location in models in which the lifetime of the LLPs is too small to reach the far location. This is, e.g., true for dark photon search if its kinetic mixing with the SM photon is not suppressed too much as discussed in [4]. On the other hand, for sufficiently large lifetime this advantage is missing, while far location suffers much less from the expected background and can additionally benefit from allowing more space for larger detector as we illustrated in the case of dark Higgs boson [5].
As we will see below, the lifetime of HNLs is typically large and they can easily overshoot detector at near location once additionally boosted. In the following, we will therefore focus exclusively on the detector placed at far location in the TI18 tunnel.
The paper is organized as follows. In Sec. II we review the basic properties of HNLs. We discuss possible production channels of sterile neutrinos at the LHC, as well as their decays that give rise to a signal in FASER in Sec. III. The main results of this study with the sensitivity reach of FASER in the searches for HNLs are presented in Sec. IV. Section V is devoted to discussion of selected scenarios going beyond the minimal seesaw mechanism that can be probed by FASER. We conclude in Sec. VI. A more detailed discussion of fragmentation functions for D and B mesons can be found in Appendix A.

II. PROPERTIES OF HEAVY NEUTRAL LEPTONS
Undoubtedly, the most important reason to add new heavy neutral leptons to the SM is to explain the neutrino oscillations [10,11] since they provide an elegant way to generate nonzero neutrino masses via the seesaw mechanism [12][13][14][15][16]. In type-I seesaw models one extends the SM by adding N neutral right-handed fermions (identified with HNLs) that couple to the SM neutrinos similarly to the coupling between left-and right-handed components of the charged leptons. Since additional right-handed fermions are SM singlets, also Majorana mass terms are allowed leading to the following Lagrangian where N I=1,...,N are the right-handed HNLs and we will assume N = 3,Φ i = ij Φ * j , Φ and L α=e,µ,τ are Higgs and lepton doublets, F αI are the Yukawa couplings between HNLs and the SM leptons and M = diag(M 1 , M 2 , M 3 ) is the Majorana mass matrix for HNLs.
After electroweak symmetry breaking (EWSB), the Higgs field gets a non-zero vacuum expectation value, v = 174 GeV, and terms in Eq. (1) proportional to Yukawa couplings generate mixing between HNLs and active neutrinos. The full 6 × 6 neutrino mass matrix reads where (M D ) αI = v F αI . One then obtains non-zero masses for the SM neutrinos after diagonalization. Similarly, the HNL mass eigenstates, N I (denoted without tilde), get small contributions from active states with mixing angles given by where typically U αI 1. HNLs inherit their interactions with the SM particles from the active neutrinos with couplings suppressed by the mixing angles. In particular the charged current interaction for the HNL is given by −(g/ √ 2) W µ¯ L,α γ µ U αI N I + h.c.. On the other hand, their masses, m N I , are mostly governed by Majorana masses M I with only small corrections from mixing. Hence, they are often referred to as sterile neutrinos or simply heavy neutrinos. If sterile neutrinos are the only particles added to the SM to explain the active neutrino masses, then the presence of at least two heavy neutrinos is required to generate two measured mass differences between the active neutrinos. At least three heavy neutrinos are required if also the lightest active neutrino wind up being massive. In this case three generation structure of the SM might also be maintained in the HNL sector.
However, not all these sterile neutrinos are necessarily within the reach of FASER. Hence, when discussing the sensitivity reach, we will focus for simplicity on a single sterile neutrino, N , with mass m N and mixing angles with the active neutrinos denoted by U eN , U µN and U τ N for the electron, muon and tau neutrino, respectively. In addition, for the purpose of presenting results, we will assume that N mixes exclusively with only one of the active neutrinos. We will therefore study the reach of FASER for three scenarios in which the mixing is dominated by only one non-zero value of U N in each case, where = e, µ, τ . On the other hand, the reach for scenarios in which more than one HNL has mass and couplings in the region of the parameter space covered by FASER, or when the mixing between N and active neutrinos has a more complicated pattern, can be deduced from our results based on the plots presenting the number of events in each of the three cases considered by us. The effective parameters that we vary are then the following: The allowed region for the mass of sterile neutrinos spans over more than fifteen orders of magnitude (see, e.g., [9] for a recent review). Theoretically the most appealing, but experimentally extremely challenging, are seesaw models with M I close to the scale of Grand Unified Theories (GUT) (see, e.g., [17]) that allow to naturally fit the experimental neutrino data with Yukawa couplings F αI ∼ 1. For much lower M I , one naively expects F αI 1 which would keep sterile neutrinos beyond the reach of low-energy experiments. However, for more than one generation of HNLs, accidental or symmetry driven cancellations between contributions to the active neutrino mass matrix coming from different N I s allow larger values of Yukawa couplings to be considered in Eq. (3) [18][19][20][21][22][23]. As a result, one can effectively study HNLs with the N I masses and mixing angles as independent parameters, while still preserving the primary motivation for such models that come from explaining the active neutrino data.
In order to explain the masses of active neutrinos via the seesaw mechanism, the mixing parameters U N need to be sufficiently large. If light enough, sterile neutrinos could then thermalize in the early Universe and distort successful predictions of the Big Bang Nucleosynthesis (BBN) by contributing to the number of relativistic degrees of freedom [24]. They could also cause an additional entropy production from their late-time decays. This can be circumvented, if HNLs decay before BBN. Combining this and current experimental bounds, one obtains an effective lower limit on the sterile neutrino mass m N 140 MeV in the range of the mixing angles that can be covered by searches for LLPs. More detailed studies [25][26][27] show that in some cases m N can be lowered down to O(10 MeV).
On the other hand, BBN bounds can be evaded for even much lighter sterile neutrino provided that its mixing angles with the active neutrinos are suppressed. This gives rise to sterile neutrino dark matter (DM) candidate with mass at ∼ keV scale that was, e.g., proposed as a possible explanation to 3.5 keV excess [28,29] in the XMM-Newton data.
One can both accommodate the observed active neutrino data, as well as keep N 1 as a DM candidate, if heavier HNLs, N 2,3 , have larger mixing angles and satisfy aforementioned lower mass limit from BBN. In particular, this is true for the renown νMSM model [23]. In such a case, while effectively stable N 1 would remain outside the reach of FASER, the displaced decays of heavier sterile neutrinos could give rise to an observable signal.
One additional motivation that lies behind low-scale seesaw models is that, similarly to high-scale models, they can also explain the baryon asymmetry of the Universe (BAU) via thermal leptogenesis [30]. The lepton asymmetry in this case is generated from CP violation in out of equilibrium production and subsequent evolution of sterile neutrinos with oscillation effects taken into account. This requires low values of Yukawa couplings F αI s. As a result, the mixing angles are typically of order U N ∼ 10 −5 − 10 −3 (see, e.g., [31] for an extensive discussion) and lie beyond the reach of current experiments. A new generation of dedicated LLP searches, including FASER, is therefore needed to study this scenario.

III. HEAVY NEUTRAL PRODUCTION AND DECAYS
Sterile neutrinos can be produced at the LHC in decays of mesons and tau leptons [6], heavy baryon decays [32], as well as in production from on-or off-shell gauge or Higgs bosons (see, e.g., [33] for a recent discussion). Given our focus on very forward going light HNLs, the main production mechanism is via meson decays. The dominant contribution comes from leptonic and semileptonic decays of pseudoscalar D and B mesons, while decays of vector mesons are sub-dominant [6]. In addition, tau decays to HNLs are taken into account for scenario in which sterile neutrino mixes exclusively with the tau neutrino. Kaon decays into sterile neutrinos are kinematically allowed only in the mass range at which strong constraints exist from previous beam-dump experiments. Hence, they will not play important role in our analysis. We do not consider charged pion decays into HNLs since they only become possible for m N m π which is disfavored by strong cosmological bounds, as discussed above. Importantly, although 2-body leptonic decays of mesons are chiral suppressed in the SM, this suppression is partially overcome if sterile neutrino mass is heavy enough since its mass, m N , replaces the mass of a charged lepton in the chiral suppression factor [34].
Contributions from baryon decays are sub-dominant in the case of charmed baryons and typically also small for bottom baryons, although in some regions of parameter space they can be as large as 15% of the total production rate [32]. We will neglect them hereafter given larger uncertainties that are expected in our signal rates.

A. HNL production in meson decays
To determine momentum and angular distribution of sterile neutrinos produced in the forward region at the LHC, we first simulate meson production in high-energy pp collisions at large pseudo-rapidities. In particular, we use Monte-Carlo event generator EPOS-LHC [35], as implemented in the CRMC simulation package [36], to simulate the kaon distributions. Distributions of heavier mesons are simulated using FONLL tool [37][38][39] in which differential cross sections for c and b quark production in high-energy pp collisions are calculated. We use the CTEQ 6.6 pdfs with m b = 4.75 GeV. Subsequent hadronization is performed with non-perturbative BCFY fragmentation functions [40] for D mesons and Kartvelishvili et al. distribution with fragmentation parameter α = 24.2 for B mesons [41,42]. A more detailed        [38,44]. The list of the decay channels into sterile neutrinos that we take into account can be found in the top panels of Figure 2 where we show relevant branching fractions [6] as a function of m N for U eN = 1 and U µN = U τ N = 0. This is then combined with the fragmentation fractions in the bottom panels of Figure 2. Below kaon threshold, leptonic K ± → e ∓ N and semileptonic K 0 S,L → π ± e ∓ N decays, which are not shown in Figure 2, play a dominant role. In the case of non-zero mixing with the muon or tau neutrino, the corresponding SM leptons are allowed in the final states when kinematically not forbidden.
For m D > m N > m K , decays of D mesons give a dominant contribution as they are produced more abundantly at the LHC than B mesons. As can be seen, the dominant production mode for HNLs is then via leptonic decays of D + s since this process is not CKM suppressed (V cs ∼ 1). For the same reason, non-negligible contribution also comes from semileptonic decays of D mesons into kaons, but these processes are more phase-space suppressed with respect to two-body decays of D + s . In the case of B mesons, which becomes important for m N > m D , the dominant contribution comes from the least CKM suppressed (V cb ∼ 0.04) channels, i.e., semileptonic decays of B ±,0 to D mesons and leptonic decays of B ± c . For m N > m B − m D , we add important contributions from leptonic decays of B ± that are more strongly CKM suppressed than decays of B ± c , but are less sensitive to a detailed modeling of fragmentation into B ± c . If sterile neutrinos mix dominantly with the tau neutrino, HNL production in decay channels of D mesons is suppressed for m N > m D −m τ . In this case, provided that m N < m τ , tau decays give a dominant contribution to HNL production, where tau leptons typically originate from D s → τ ν τ decays.
The kinematic distributions for HNLs produced in B ± (D ± s ) meson decays are shown in the top central (bottom central) panels of Figure 1 for m N = 3 GeV and U eN = 2 × 10 −3 (m N = 1 GeV and U eN = 4 × 10 −4 ) while other mixing angles are set to zero. The choice of the mixing angles for both points is dictated by the current bounds that exclude values above a few ×10 −3 or a few ×10 −4 for m N above or below the D meson threshold, respectively. For integrated luminosity of 3 ab −1 at the 13 TeV LHC one can obtain about 10 6 (10 8 ) sterile neutrinos with these masses and mixing angles.

B. HNL decays in FASER
Once produced in the forward direction, HNLs can travel long distances before decaying as their decay width is suppressed by the square of the mixing angle, |U N | 2 . In the left panel of Figure 3 we show typical decay lengths cτ N , where τ N is the sterile neutrino lifetime, assuming non-zero mixing with the electron neutrino with |U eN | 2 = 10 −7 . The lifetime in case of mixing with the muon or tau neutrino is of similar order with somewhat more longlived HNLs obtained in the latter scenario due to kinematical suppression of modes with the tau lepton in the final state. Importantly, as shown in the central panels of Figure 1, HNLs produced at the LHC in the forward direction are typically boosted what further increases their decay length in the lab frame. As can be seen in Figure 3, it often exceeds 1 km with lower values possible only for m N m D . Hence, only a small fraction of produced HNLs effectively decay at the location of FASER. However, as we will see below, even a small detector might find enough number of HNL decays to probe interesting part of the parameter space of this model.
The probability of decay within the detector volume is given by whered = γ β c τ N , θ N is the angle between the sterile neutrino momentum and the beam axis in lab frame, R is detector radius, L min corresponds to the position of the detector and L max = L min + ∆ where ∆ is the length of FASER along the beam axis. Since typicallȳ d L max , the probability of decay within the detector volume scales linearly with ∆ and is  inversely proportional to the decay lengthd as shown in Eq. (5). Simple Eq. (5) is modified in case the particles decaying into HNLs have non-negligible lifetimes and, therefore, can travel sizable distances before decaying. In particular, it is true for K 0 L and K ± decays into HNLs. The former is required to decay before hitting the first neutral absorber, while the latter before reaching first magnets deflecting their trajectories.
We consider two representative detector setups at the far location: L max = 480 m, ∆ = 10 m, R = 20 cm, 1 m.
A small detector with R = 20 cm is technologically more appealing and is large enough to achieve the full reach of FASER in the case of search for dark photon [4] and axion-like particles [46]. On the other hand, larger detector with R = 1 m is desired if new physics particles are produced dominantly in heavy meson decays like, e.g., dark Higgs bosons [5]. In the right panels of Figure 1 we show the kinematical distribution of sterile neutrinos that decay within the volume of FASER for the same points in the model parameter space as mentioned above. As can be seen in the top panel, in case of m N = 3 GeV only high energy HNLs with E N few hundred GeV can reach the detector. On the other hand, as illustrated in the bottom panel, HNLs with lower masses and couplings have larger lifetime and therefore can survive until they reach FASER even if they have energies below 100 GeV. An approximate region of the parameter space in the (m N , U eN ) plane that can be probed by FASER can be deduced from the right panel of Figure 3 where we show the contours of constant decay length of 1 TeV sterile neutrinos. We can see that most parts of the parameter space correspond to the long lifetime limitd L max in which the HNL tends to overshoot the detector and the event rate is limited by the detector length. On the other hand, at large masses m N > m D and mixings U eN ∼ 10 −2 the decay length becomes short d < L max and only highly boosted particles can reach FASER.
In Figure 4 we show the main decay channels of sterile neutrinos for each of the three scenarios with non-zero mixing angles that we consider. The leptonic channels include various ν l + l − final states, as well as invisible decay modes into three SM neutrinos. For m N > m π , various 2-body decays into a pair of a lepton and a meson are possible which we implement following [6,7]. In particular, the green lines for m N 1.2 GeV correspond to combined branching fraction to various light unflavored mesons including π 0,± , η, ρ, η , ω and φ, as well as to strange mesons, K ( * ),0,± . We further assume that for m N 2 GeV one can consider quarks as degrees of freedom in the final state in decays N → (ν/l)q 1q2 [31]. For intermediate values of the HNL mass, 1.2 GeV m N 2 GeV, we interpolate hadronic decay widths requiring that branching fractions for all decay modes to sum up to unity. The actual values of all branching fractions in this mass regime should be treated less robustly which is denoted by the use of dashed lines.

IV. SENSITIVITY OF FASER TO HNLs
The expected signal from sterile neutrinos decaying in FASER consists of two simultaneous high-energy charged tracks, which originate from a vertex inside the detector and whose combined momentum points into the direction of the IP. While detailed background analysis goes beyond the scope of this paper, one can argue that natural and infrastructure-based shielding, as well as specific properties of the signal, should allow to disentangle it from background [4]. In particular, SM particles produced abundantly at the IP would be either deflected by the LHC magnets or stopped by neutral absorbers before they reach FASER. On the other hand, energetic particles produced in beam-gas collisions in the beam pipe close to the detector would effectively loose their energy in the rock and concrete separating the side tunnel from the main tunnel. The kinematic features of the signal, such as the directionality of the charged particles observed in FASER, would provide an additional handle that could help to discriminate between these backgrounds and the signal events, as well as to reject events with cosmic-ray origin. In the following, we will for simplicity assume zero background when discussing FASER's sensitivity.
FASER's sensitivity reach can also be affected by the finite detection efficiency. In particular, we correct the expected number of events to take into account sterile neutrino decays into three SM neutrinos. In the following, we will assume 100% efficiency for other decay channels, while a more detailed analysis of this effect is postponed to a later dedicated study.
In Figure 5 we show the sensitivity reach of FASER to search for sterile neutrinos that mix only with the electron neutrino, where the gray bands correspond to the current exclusion limits. In the left panel we show the sensitivity of R = 1 m detector for scenario with U eN = 0 decomposed into contributions from kaon (blue), D meson (green) and B meson (red) decays. Typically lighter meson decays dominate when kinematically allowed since they are produced more abundantly at the LHC. The solid lines correspond to a fixed number of signal events, n = 3, 10, 10 2 , 10 3 , 10 4 . As can be seen, one can expect up to ∼ 1000 events with HNL origin in FASER for m N m D , while ∼ 10 signal events for sterile neutrinos produced in D meson decays. Clearly, possible changes in the number of signal events by a factor O(1) from a finite efficiency of the detector or non-negligible background would have a mild impact on the reach especially in the case of B meson decays.
In the right panel of Figure 5 we compare the sensitivity of R = 20 cm and R = 1 m detector for 3 ab −1 integrated luminosity. As can be seen, the reach for R = 1 m is significantly improved with respect to the one obtained for a smaller detector. This is not a surprise since HNLs, similarly to dark Higgs bosons [5], are typically produced in heavy meson decays and, therefore, they are less collimated around the beam axis than, e.g., dark photons produced in pion decays [4]. The reach could be further improved, especially for m N m D , for even larger radius, as can be deduced from the bottom right panel of Figure 1. Similarly, the case of mixing with the muon neutrino and tau neutrino are shown in Figure 6 and Figure 7, respectively. As can be seen, FASER will probe parts of the yet unconstrained region of the parameter space in each of these scenarios. In the electron and muon case, current bounds below kaon threshold are stronger than the reach of FASER. However, FASER capability is much improved once m N grows and the typical decay length of boosted sterile neutrinos becomes comparable to the distance between the IP and the detector (compare Figure 3). On the other hand, the scenario with non-zero mixing with the tau neutrino is much less constrained. This leads to even larger region of the parameter space for HNLs that has not been probed yet and is characterized by an excellent discovery potential in FASER with possibly even tens of thousands of expected signal events for 3 ab −1 integrated luminosity.
In the range of m N relevant for FASER, the most stringent bounds on sterile neutrino mixing with the electron neutrino or muon neutrino come from past and present beam-dump experiments at CERN (PS191 [52], CHARM [53] and NA62 [54]) and IHEP-JINR [55]. For m N 2 GeV the strongest limits come from search for B decays into HNLs at Belle II [56] and from limits on the Z boson decays into HNLs from the LEP data collected by the DELPHI collaboration [57]. For mixing with the muon neutrino other important bounds come from search for a double-peak structure in K → µν decays [58] and the NuTeV beam-dump experiment [59]. In the case of mixing with the tau neutrino, current limits are much weaker with the leading bounds coming from CHARM [60] and DELPHI [57] collaborations. If the mixing angles become low enough, strong bounds from BBN [25][26][27] constrain the parameter space of HNLs from below. In scenario with non-zero mixing with the electron neutrino, U eN = 0, other important bounds come from null searches of the neutrinoless double beta decay, denoted as 0νββ.
The most stringent current limit on the 0νββ decay half-life T 0ν 1/2 1.07 × 10 26 yr comes from a combined analysis of the Phase-I and Phase-II data acquired by the KamLAND-Zen experiment [61]. This can be translated (see, e.g., [62,63]) into an approximate limit on the mixing angle |U eN | 2 /m N 2.1 × 10 −8 GeV −1 . However, as discussed above, more than one sterile neutrino is required for the seesaw mechanism to generate correct active neutrino masses. In addition, some of the sterile neutrinos can be lighter than the typical momentum transfer q ∼ 100 MeV in 0νββ process provided they interact weakly enough so as not to alter BBN. In this case, additional cancellations between various contributions to the 0νββ rate may occur and effectively weaken the corresponding bound [64]. It is therefore not shown in Figure 5.
For comparison, in Figure 5 we also show the expected sensitivities for sterile neutrino searches in the proposed SHiP detector [47], the DUNE [48] and NA62 [49] experiments and in the LHC searches for prompt lepton plus a single displaced lepton jet [50]. As can be seen, FASER has comparable sensitivity to NA62 for sterile neutrinos produced in D meson decays. This region in the parameter space for both electron and muon mixing will also be probed by the future DUNE facility for which we show results following [47]. They correspond to the analysis performed for the 5 years of data-taking by the 30 m long LBNE near detector with 5 × 10 21 protons on target and assuming a normal hierarchy of neutrinos.
Above the D meson threshold, the abundant production of forward B mesons at highenergy pp collisions at the LHC works in favor of FASER. This allows FASER to probe the region between the aforementioned planned searches and the expected reach of the LHC searches [50] as shown for scenario with U µN = 0 in the bottom left panel of Figure 5. The LHC searches will also be sensitive to other scenarios with non-zero U eN or U τ N . However, as discussed in [50], dedicated analysis of relevant backgrounds and tagging efficiency would have to be performed by experimental collaborations to obtain solid results. Therefore, they are not shown in Figure 7.
In case of non-zero mixing with the tau neutrino, U τ N = 0, HNL production in D meson decays is restricted to m N < m D − m τ . The main contribution comes from the decay D + s → τ N . Larger masses m D − m τ < m N m τ can be probed by hadronic and leptonic tau decays, τ → N π, N K + , N ρ + and τ → N ν . Here the τ s are mainly produced via the decay D + s → τ ν. This region of the parameter space can also be partially probed by searches based on τ production in B factories [51] and their subsequent decays into sterile neutrinos as shown as magenta line Figure 7.
The strongest expected reach in the range of the HNL masses probed by FASER could be obtained by the proposed SHiP detector. In order to compare the reach of FASER to the one of SHiP, it is useful to note that to a good approximation the number of expected events scales like where in the second step we used Eq. (5) and for simplicity we neglected the geometrical acceptance of the detector. The expected number of mesons produced for a given experiment is denoted by N mes , while E mes is the typical energy of mesons that give rise to sterile neutrinos decaying within the volume of each detector. For both considered experiments the relevant numbers read Combining Eq. (7) and Eq. (8) one obtains for SHiP only about ∼ 5 times larger number of events with B-meson origin. Given that in the regime of large lifetime the number of expected events scales like N ev ∼ |U N | 4 , one obtains only slightly better sensitivity (by a factor of 2) of SHiP to sterile neutrinos produced in B meson decays in the (m N , U N ) plane. In addition, for somewhat heavier B ± c mesons, FASER can gain even more in meson production from much larger collision energy than SHiP. As a result, the sensitivity of both detectors is comparable at the largest accessible values of m N ∼ 4 GeV. The loss of sensitivity for this large masses and increasing mixing angles to values U N 10 −3 can be understood since then the sterile neutrino lifetime becomes too low for HNLs to reach FASER.
On the other hand, production of D mesons is less suppressed for SHiP center-of-mass energy 27 GeV and therefore larger luminosity allows to obtain ∼ 10 3 −10 4 more HNL events with D-meson origin in SHiP than in FASER. Hence, the reach in the mixing angle is worse for FASER by about an order magnitude as can be seen in Figure 5. This is also true for sterile neutrinos produced in τ decays, as seen in Figure 7, since they typically come from D meson decays.

V. BEYOND MINIMAL SEESAW
In previous sections we have analyzed in details FASER's sensitivity reach to HNLs described effectively by only their mass and the mixing angles with the active neutrinos. We will now discuss how this search can be relevant to selected more complex scenarios going beyond the minimal seesaw mechanism with only right-handed neutrinos added to the SM. In particular, the reach for historically very important Left-Right symmetric models is discussed in Sec. V A. In Sec. V B we illustrate possible connection of the sterile neutrino search in FASER to DM and the baryon asymmetry of the Universe.

A. LR symmetric models
In the Left-Right (LR) symmetric models [65][66][67] the SM gauge group is extended to SU (2) L × SU (2) R × U (1) B−L with gauge couplings of left and right SU (2) groups denoted by g L and g R , respectively. As a consequence of the so-called left-right symmetry, additional right-handed doublets of quarks and leptons that appear in the model are constructed analogously to the left-handed doublets of SU (2) L . In particular, this naturally implies the existence of three right-handed neutrinos that complete the right-handed lepton doublets and can lead to non-zero masses of the active neutrinos via the seesaw mechanism. The LR symmetry is broken to the SM gauge group, SU (2) L × U (1) Y , at some energy above the scale of EWSB. This leads to additional charged and neutral gauge bosons denoted by W R and Z R , respectively. Current exclusion bounds imply m W R m W L , while Z R is even heavier and we will henceforth neglect its impact.
Right-handed neutrinos in this model can be produced and subsequently decay in various processes involving W R exchange. Assuming, for simplicity, that g L = g R and given that mixing matrix in the right-handed quark sector is expected to follow closely the left-handed one [68], the matrix elements of the relevant processes are rescaled with respect to the SM neutrinos by (m W L /m W R ) 2 . Hence, phenomenology of the right-handed neutrinos in LR symmetric models closely resembles the above discussion of HNLs with a simple substitution  4 . This allows to easily translate limits on HNLs from various intensity frontier searches into the limits in the (m W R , m N ) plane (see, e.g., [69]). Importantly, sterile neutrinos in the LR symmetric models with sizable mixing with the electron neutrino are strongly constrained by null searches of neutrinoless double beta decay [70,71]. For this reason we present the results for LR sterile neutrino from muon doublet assuming that other mixing angles are suppressed. In Figure 8 we present sensitivity reach of FASER in the (m W R , m N ) plane for such scenario. As can be seen, FASER could provide limits reaching m W R 7 TeV for low mass of N that would be complementary to the region of the parameter space probed in search for W R → µN → µµjj events in the high-p T searches at the LHC. In particular, we show current such limits obtained by the CMS collaboration for √ s = 13 TeV and 35.9 fb −1 integrated luminosity [72], as well as by the ATLAS collaboration for √ s = 8 TeV and 20.3 fb −1 integrated luminosity [73]. Note that the ATLAS search also considers a topology with one jet in the final state extending their reach for light and highly boosted m N . Additionally, further complementarity is achieved once searches for multitrack displaced vertices at the LHC [74,75] are taken into account that are sensitive to scenarios with m N 10 GeV. We also show for comparison the projected sensitivity of the proposed SHiP detector obtained based on [47] which could further extend the future bounds up to m W R ∼ 20 TeV for a narrow range of the HNL masses. Other bounds, which are not shown in Figure 8, from searches for HNLs that we have discussed in Sec. IV can also constrain the parameter space of the LR symmetric models for low values of m W R . B. Sterile neutrinos from the mirror sector FASER can also be sensitive to scenarios employing other proposed mechanisms of generating the masses of the active neutrinos that go beyond type-I seesaw. Such possibility with additional connections to DM and leptogenesis has been discussed in [76] for the SM extended by its mirror sector. Notably, models with the mirror sector of the SM have recently received much attention also in the context of studies of neutral naturalness as a solution to the hierarchy problem. In particular, a prototypical such scenario is the Twin Higgs model [77] with the SM duplicated in the mirror SM-singlet sector. Needless to mention that neutral naturalness becomes increasingly motivated given lack of discovery of any colored BSM particles at the LHC.
In model considered in [76] it is further assumed that both the SM and the mirror sector are connected via set of heavy right-handed Majorana neutrinos, N I , that have Yukawa couplings to the SM leptons, as well as to their mirror counterparts where prime corresponds to the mirror sector. We will denote the active neutrino states in the SM and the mirror sector by ν and ν , respectively. In addition, for consistency, the set of Higgs bosons of the model has been extended to include Higgs doublets, H u,d (H u,d ), and triplets, ∆ (∆ ), for the SM (mirror) sector. The latter couple to lepton pairs via Y ∆L c ∆ L and Y ∆L c ∆ L . The neutral components of both triplets can get non-zero vev which gives raise to Majorana (type-II seesaw) contributions to the neutrino mass matrix for ν and ν that are denoted by µ and µ , respectively.
Leptogenesis in both sectors can be driven by out-of-equilibrium decays of the lightest of the right-handed neutrinos or by resonant processes if two right-handed neutrinos are mass degenerate. It is then transferred to the baryon asymmetry by sphaleron transitions. In particular, baryons of the mirror sector can play a role of DM. Correct ratio between baryon and DM abundances is set by the ratio between the SM and mirror proton masses. The latter is driven to the value ∼ 5 m p by assuming large scale of the EWSB in the mirror sector, v mir /v SM ∼ 10 3 − 10 4 that makes mirror sector u and d quarks heavy.
After diagonalizing the neutrino mass matrix, the SM active neutrino masses have both type-II seesaw contribution µ and the inverse seesaw contribution that can be written as µ (v SM /v mir ) 2 . 1 On the other hand, ν masses in the mirror sector are mostly driven by type-I seesaw mechanism. Typically two of them are in the ∼ GeV scale for the parameters of the model that allow successful leptogenesis and correct DM abundance [76,79]. In addition, in the mass eigenstates, non-zero mixing between ν and ν appears of order v SM /v mir ∼ 10 −4 − 10 −3 . As a result, ν s play a role of the light sterile neutrinos with properties resembling those of HNLs discussed above. Importantly, their typical masses and mixing angles lie in the region of the parameter space covered by FASER which should, therefore, be able to test this interesting cosmological scenario.

VI. CONCLUSIONS
Searches for light long-lived particles have recently become widely accepted as one of the most important aims to achieve in the next generation of dedicated experiments looking for BSM physics. Among many models of new physics that can be probed in such searches, one of the most extensively studied scenarios predict heavy neutral leptons with masses of order GeV and small mixing angles with the active neutrinos. In particular, HNLs appear naturally in the context of the seesaw mechanism that generates non-zero neutrino masses in the SM. Although typically seesaw scenario predicts sterile neutrinos with large masses of order the GUT scale, it is also possible to obtain light HNLs in theoretically appealing models.
In this study we have analyzed the expected sensitivity reach to HNLs of the newly proposed detector to be placed along the LHC beam axis, named ForwArd Search ExpeRiment, or FASER [4,5]. This analysis is complementary to our previous studies for dark photons that are mainly produced in pion decays and for dark Higgs bosons that typically come from B meson decays. In the case of HNLs also decays of D mesons play an important role, while B mesons start to dominate only when other channels are kinematically forbidden.
We compare the sensitivity reach of FASER in the HNL parameter space to other proposed detectors. In particular, we show that for sterile neutrinos produced mainly in B meson decays, even a relatively small cylindrical detector with ∼ 10 m length and R 1 m radius can have comparable reach to the proposed SHiP detector, while other experiments typically loose their sensitivity in this range of the HNL mass. In this mass regime the current bounds in the mixing angles could be improved by about an order of magnitude. On the other hand, the reach in the mixing angle for HNLs produced in D meson decays is comparable to the one expected from the NA62 experiment and can be surpassed by the SHiP and DUNE due to their much larger luminosity.
However, even in the case of D meson decays FASER can probe interesting region in the HNL parameter space. In particlar, as discussed in Sec. V B, it corresponds to the attractive cosmological scenario in which sterile neutrinos emerge from the mirror sector of the SM with heavy right-handed neutrinos as mediators. In this scenario, search for sterile neutrinos in FASER could shed light on other important questions in contemporary physics including the nature of DM and the origin of baryon asymmetry of the universe. It could also be related to increasingly popular solution to the hierarchy problem via neutral naturalness.
Importantly, looking for a signal from light sterile neutrino decays in FASER will play a complementary role to searches for heavier HNLs in high-p T oriented experiments at the LHC. We illustrate this in Sec. V A for the case of the left-right symmetric models that provide a natural extension of the SM gauge group to explain the asymmetry between the left and right fermions in the SM.
Last, but not least, in scenarios with sterile neutrinos emerging from the mirror sector of the SM, the photon of the mirror sector mixes with the SM photon [76,79]. It can then effectively play a role of the dark photon. Its mass and mixing parameter would then typically be within the reach of FASER discussed in our previous study [4]. In this case, one expects to simultaneously see events with the HNL and dark photon origin with possibly distinct final states of their decays. As discussed in Sec. III, the dominant contribution to production of HNLs in the context of FASER comes decays of pseudoscalar D and B mesons. To determine the D meson fragmentation fractions and functions we follow [80]. The fragmentation tree for charm quarks is shown in Figure 9.
The charm quark fragments into mesons with a probability of P m = 0.89 or baryons with P b = 0.11. The baryonic mode is dominated by the production of Λ c with an 14% contribution from charmed-strange mesons Ξ c and Ω c [81]. The meson are either produced in the pseudo-scalar or vector state with production fractions P P = 0.39 and P V = 0.61 respectively [80]. The fragmentation fraction into the strange meson states D ( * )+ s is given by, p s = γ s /(2 + γ s ) ≈ 0.12, where we used strangeness suppression factor γ s = 0.259 [80]. For u and d quarks, the fragmentation fractions are obtained from p u,d ≈ (1 − p s )/2 ≈ 0.44 provided that the ratio of neutral and charged mesons is equal to R u/d = 1.02 ≈ 1 [80]. The vector meson states will then decay into the pseudo-scalar state and either a photon or pion. The dominant decay modes of the D * + d meson into D 0 or D + d have branching fractions B(D * + d → D 0 π + ) = 0.68 and B(D * + d → D + d (π 0 /γ)) = 0.32, respectively [82]. These numbers are then combined into pseudo-scalar and vector fragmentation fractions, f P and f V , for the pseudo-scalar mesons D 0 , D + d and D + s as shown on the right side of Figure 9. The pseudo scalar and vector fragmentation fractions f P,V for each meson have different non-perturbative fragmentation functions F P (z) and F V (z), respectively, which depend on the momentum fraction z. We employ the BCFY fragmentation functions [40] with default FONLL settings [38]: m c = 1.5 GeV, N = 5 and r = 0.1. The total fragmentation function for each meson is given by where D i = D 0 , D + d and D + s are the pseudo-scalar D mesons for which we obtain the kinematic distributions.

B mesons Fragmentation
To simulate B-meson fragmentation we use a non-perturbative fragmentation function which follows distribution by Kartvelishvili et al. [41,42] with fragmentation parameter α = 24.2. At this stage we do not differentiate between various B-meson final states.
To determine the number of various B hadrons we employ fragmentation fractions f u , f d , f s , f c and f Λ that correspond to B + u , B 0 d , B 0 s , B + c and Λ b hadrons in the final state. For simplicity, we include all baryonic states into f Λ so that above fragmentation fractions sum up to unity. We then employ the LHCb results for f i [83,84] that correspond to continuum bb production mode Note that we only take into account results from LHCb, which use the continuumbb production mode. These results agree very well with the fragmentation spectrum measured at Tevatron. The fragmentation fractions measured at BaBar and Belle at the Υ-resonance and LEP at the Z-resonance differ slightly due to the resonant production mode. We can solve for the individual fragmentation fractions f i and obtain