Discovery potential of FASER$\nu$ with contained vertex and through-going events

FASER$\nu$ is a newly proposed detector whose main mission is to detect neutrino flux from the collision of the proton beams at the ATLAS Interaction Point (IP) during the run III of the LHC in 2021-2023. We show that this detector can also test certain beyond standard model scenarios, especially the ones in which the neutrino interaction with matter fields can produce new unstable particles decaying back into charged leptons. Models of this kind are motivated by the MiniBooNE anomaly. We show that, if the new physics involves multi-muon production by neutrinos scattering off matter fields, including the neutrino flux interactions in the rock before the detector in the analysis (i.e., accounting for the through-going muon pairs) can significantly increase the effective detector mass and its sensitivity to new physics. We propose a concrete model that can give rise to such a multi-muon signal.


I. INTRODUCTION
Neutrinos are copiously produced at the interaction points of the LHC experiments. However, their detection at the main detectors of the LHC (i.e., CMS, ATLAS, ALICE and LHCb) is not possible because of the large background from other particles produced at the Interaction Point (IP). To detect high energy neutrino flux from the IP, the FASERν experiment [1] has been proposed. It will take data during run III of the LHC at a distance of 480 m upstream of the ATLAS interaction point. The FASERν detector is composed of layers of Tungsten interleaved with emulsion film with a total size of 25 cm×25 cm×1.3 m and a total mass of 1.2 ton. As the neutrinos travel in a straight line and traverse 10 meters of concrete and 90 meters of soil, they can reach the detector.
As shown in [1][2][3][4], FASERν can also search for new physics. We focus in particular on the muon channel. FASERν can detect charged lepton produced in the Charged Current (CC) interaction of all three flavors of neutrinos inside the detector. In principle, the muons produced by ν µ CC interactions in the rock and concrete before FASERν can also give rise to a signal. However, the background of muons from IP makes using these muons challenging.
For the first time, in this paper, we show that in the presence of multi-muon events due to new physics the effective mass of the detector can be increased by a factor of 4 by including events starting in the rock or concrete before the detector. This effect is similar to the enhancement of the effective volume of ICECUBE to detect ν µ by searching for the through-going muons from below. We point out that, if FASERν is equipped with a scintillator detector that can record the arrival time of the muons comprising a multi-muon event, the background from the pile-up of the muons penetrating from IP as well as from the charged current interaction of ν µ in the rock can be made negligible.
In this paper, we first show that with this technique it is possible to considerably enlarge the effective mass of the detector by registering the multi-muon events produced outside the detector. We then introduce a model that leads to a signal of charged lepton pairs. The model, which is also motivated by the MiniBooNE anomaly [5,6], includes a new right-handed neutrino, N , with a mass in the range of few 100 MeV to ∼ 10 GeV and a dark photon, Z ′ with a coupling converting ν and N to each other. The well-known technique to obtain such off-diagonal Z ′ interaction is to mix the new fermions with the Standard Model (SM) neutrinos. Since there are strong upper bounds on the mixing, obtaining large couplings for the off-diagonal Z ′ interaction is usually a challenge but in the appendix, we introduce a new trick to circumvent the bounds and obtain a coupling as large as 0.01.
With this coupling, The scattering of the neutrino flux off nuclei will create N . If N is heavier than Z ′ , it can decay into Z ′ which in turn decays into e − e + and, for heavy enough Z ′ , also into µ − µ + . If N is lighter than Z ′ , it will go through three-body decays, producing e − e + or µ + µ − pairs. By studying the energy-momentum and tracks of the final charged leptons originated from a vertex inside the FASERν detector, further information about the lifetime and masses of the new particles can be derived. We also show how we can build a variation of the model in which the N decay leads to a four muon signal instead of a dilepton one. We demonstrate how including the through-going four-muon signals can extend the capability of FASERν to probe smaller values of the neutrino coupling to the new particles.
This paper is organized as follows. In sect. II, we show under what conditions including the neutrino interactions in the rock and concrete can enlarge the effective volume of the detector. In sect. III, we outline the basics of the model that can lead to neutrino upscattering to N and its subsequent decay that can yield a multilepton signal. The details of the model, with emphasis on the trick to obtain large upscattering cross-section for neutrinos and on obtaining the µμµμ signal, are presented in the appendix. In sect. IV, we discuss the decay modes of N and their signature at FASERν. In sect. V, we show how to calculate the production rate of N from the ν interaction and how to compute the number of multi-lepton events that can be registered at FASERν. The results on the number of events and the bounds that can be obtained by FASERν and its upgrades are presented in sect. VI. A summary and a discussion on the extension of the results to other theoretical models are given in sect. VII with an emphasis on the role of time recording scintillator.

II. THROUGH-GOING AND CONTAINED VERTEX MULTI-MUON EVENTS AND THEIR BACKGROUNDS
In this section, after reviewing backgrounds for multimuon events and suggesting strategies for reducing them, we discuss under what conditions the through-going multi-muon events can be invoked to discover new physics. The only particle other than neutrinos that can travel in the rock up to the detector is the muon. Muons will lose energy in the rock so reconstructing their energy-momentum at the detector will not give much information. However, their directions remain fairly faithful to their initial direction at the production. In fact, at these energies and for these distances the deviation from the original direction will be less than 0.5 • (see the last paragraph of page 5 of Ref. [7]). The angular resolution of the detector is much better than this value [1] and at the level of sub-milliradians so we can practically neglect the uncertainty due to the angular resolution. Consider a muon produced at a distance of d = 14 m from the detector. The two lines connecting the production point to the center and the edge of the detector of size 25 cm × 25 cm subtend an angle of [((0.25/2)/14) × (180/π)] • = 0.5 • . Thus, roughly speaking, for the muons produced in the vicinity of 14 m from the detector (despite the deviation of the direction of the muon during the propagation) we can distinguish that they have originated in a distance between (7 m to 14 m) from the detector along the beam and are not coming directly from the Interaction Point.
According to [1], during run III of the LHC (2021-2023), O(10 9 ) muons will reach FASERν from the IP or the interaction of the neutrino flux in the rock before the detector. The probability that two separate events out of these 10 9 muons pile up in a single bunch crossing (∼ 1 ns) is negligible. Thus, recording the timing of the multi-muon event will significantly reduce the background. FASERν is an emulsion detector that cannot record time but if, as it is suggested in [1], the detector is equipped with a scintillator plate, recording the timing of the arrival of the muons with nanosecond accuracy will become a reality. Still the dimuon events can originate from pair production from photons [8] or from charm production via ν µ CC interaction (i.e., ν µ + nucleon → µ + c + X) and the subsequent leptonic decay of the charm. The number of such charm-induced dimuon events originating outside the detector can exceed B =1000 events during run III of the LHC. As we shall discuss in sect. IV, the background events with a vertex inside the detector can be suppressed by the measurement of the energy-momentum of the dimuon or by considering the event topology but these methods cannot be used for suppressing the background for dimuon events originating in the rock. FASERν can declare discovery of new physics with dimuon signal from the rock only if the number of events is much larger than the root of the number of the background events; i.e., the number of dimuons from new physics is much larger than √ B = 30.
If the new physics signal instead of dimuon is three-muon or four-muons, we do not need to worry about the background. If FASERν records more than two muon events within 1 ns, we can make sure that they originate from new physics even if we do not reconstruct the directions of the muons. For such signal events, the effective mass of the detector will be enhanced by including the events originating in the rock before the detector.
In the next section and the appendix, we shall develop a model in which one or two pairs of forward going muons are produced by the interaction of the neutrino flux in the rock and concrete before the detector. We then demonstrate how invoking the through-going muon pairs will increase the sensitivity of FASERν to search for new physics.

III. THE MODEL
Refs. [5,6,9] propose a novel idea to explain the low energy e − excess observed in the MiniBooNE experiment. It is based on introducing a sterile neutrino of mass of a few 100 MeV mixed with ν µ and a new light U (1) gauge boson, Z ′ , coupled both to the quarks and to the neutrinos with a coupling of form Through this coupling, ν µ (orν µ ) can scatter off nuclei converting to N (or toN ). In the models introduced in Refs. [5,6], the coupling of Z ′ to the quarks emerges because of the kinetic mixing between Z ′ and the photon so Z ′ also couples to e − e + . As a result, the produced N subsequently decays producing an electron-positron pair via N → ν(Z ′ ) * → νe − e + 1 As shown in [5], in the MiniBooNE experiment, these pairs can be misidentified as charged current interactions of ν e orν e . In the case that m Z ′ > m N , as is assumed in [5], the N decay will be a three-body decay mediated by an off-shell Z ′ . On the other hand, in the case m Z ′ < m N , as is assumed in [6], N can decay into Z ′ ν and Z ′ subsequently decays into the e − e + pair.
In the model presented in [5], N is charged under the new U (1) gauge symmetry so we would have a coupling of g ′N γ α N Z ′ α where g ′ is the gauge coupling which can be of order of 1. A mixing between N and ν µ leads to the coupling of form (1) with g νN = g ′ U µ4 . To be precise ν µ in Eq. (1) is not a flavor (or electroweak) eigenstate but a linear combination of light eigenstates as Since in this model g νN is proportional to the mixing (g νN = g ′ U µ4 ), for a given m N , the bound on the mixing will be translated into a bound on g νN .
Let us enumerate the bounds on the mixing of a sterile neutrino, ν 4 of a mass of m 4 with active neutrinos. For m 4 < m π − m µ ≃ 30 MeV, the search for π → µν 4 in the PIENU experiment [10] sets an upper bound of 10 −3 on U µ4 . On the other hand, for 70 MeV < m 4 < m K − m µ the searches for K → µν 4 by KEK [11], by E949 [12] and by NA62 [13] yield a very stringent bound on g νN . For heavier ν 4 , the bound comes from NuTeV [14]. All the bounds are summarized in Ref. [15]. However, for 30 MeV < m 4 < 70 MeV, the strongest bound on U µ4 , which comes from studying the muon decay spectrum [16], is relatively relaxed and is of order of few×10 −2 [15].
As we demonstrate in the appendix, it is possible to build models for the coupling of form in Eq (1) with nonzero g νN even with a vanishing mixing between N and ν µ . Throughout the paper, we, therefore, take g νN to be independent of U µ4 and study its signatures at FASERν. In addition to (or instead of) ν µ , ν e and/or ν τ could have a coupling of form (1) to Z ′ and N . We however focus on the case that the coupling is exclusively to ν µ because of two reasons: (i) The flux of ν µ at FASERν is higher than those of ν e and ν τ ; (ii) To keep connection with the solution to the MiniBooNE excess. However, similar analysis and discussion can be repeated for the ν τ or ν e coupling to N and Z ′ .
The basis of the toy model described in the appendix is introducing two sterile neutrinos N ′ and N with an offdiagonal coupling of the form g ′N γ µ N ′ Z ′ µ . The mixing of N with active neutrinos is taken to be zero or negligible. N ′ has a mass in the range (30 MeV, 70 MeV) so it can have a relatively large mixing with ν µ leading to a relatively large g νN . We take N to be heavier than N ′ so that N can decay into N ′ . If the coupling of Z ′ to the SM fermions is through its mixing to the photon, ǫ, it will not couple to νν and as a result, at the tree level, the decay modes N → Z ′( * ) N ′ → N ′ νν are closed but we can have As discussed in the appendix, Z ′ may also couple to a pair of new scalars, a andā, whose decays produce a pair of µμ. In this case, the N decay can produce the background-free signal of two pairs of µμ N → N ′ Z ′( * ) → N ′ aā followed by a → µμ andā → µμ. (2) The N ′ particle can decay via the electroweak interactions as 6 which is much larger than the distance between FASERν and the interaction point.
We now have a model for the coupling of form in Eq. (1). In order to produce N via scattering of ν off the nuclei, we also require Z ′ to couple to the quarks. Let us take the coupling of Z ′ to the standard model fermions, f as As discussed in [5], this can be achieved by the following kinetic mixing between the photon field strength F µν and the field strength of the new U (1) gauge symmetry, Going to the canonical basis in which the kinetic terms are diagonal and properly normalized, all the charged fermions of the SM, f , will obtain a coupling to Z ′ proportional to their electric charge: For 20 MeV < m Z ′ < 10 GeV, the strongest bound on ǫ comes from BABAR with ǫ < 7 × 10 −4 [17] (see also [18,19]). For such values of ǫ, the Z ′ decay in the early universe will take place long before the neutrino decoupling so that there would be no deviation from the standard big bang nucleosynthesis prediction. If Z ′ is lighter than a few hundred MeV, it can be produced via γ + e − → e − + Z ′ inside the supernova core and become thermalized until the temperature drops below the Z ′ mass. Since the Z ′ mean free path, as well as its decay length, is much shorter than the size of the core, the Z ′ production will not dramatically affect the supernova core evolution. If N is heavier than a few 100 MeV, it cannot be produced inside the supernova core but it can change the mean free path of ν µ andν µ in the inner core via ν µ +ν µ → Z ′ + Z ′ . However, at the neutrinosphere with a temperature lower than the Z ′ mass, this process cannot take place so the neutrino emission duration will not be significantly prolonged. N ′ can be produced via ν µ Z ′ → N ′ Z ′ with a rate of ∼ (g ′ g νN ) 2 T 5 /(4πm 4 N ) and via the weak interactions with a rate of ∼ G 2 F T 5 |U µ4 | 2 /(4π). The produced N ′ will be thermalized by scattering Z ′ N ′ → N ′ Z ′ via a t-channel N exchange with a mean free path in the inner core smaller than that of the active neutrinos so it cannot transfer the energy of the core outside or affect the convection in the core. Once the core is depleted of neutrinos and the temperature drops below the Z ′ mass, N ′ cannot be reproduced and it decays within ∼ 0.01 sec emitting ν µ . With the current precisions, the supernova constraints cannot rule out the scenario but it can alter the ν µ spectrum which in the future can be tested. Unless stated otherwise, throughout this paper, we shall assume that the coupling of the Z ′ to the quarks and charged leptons is through kinetic mixing between Z ′ and the photon so the relation in Eq. (3) holds.
Another option that we shall entertain in this paper is gauging the anomaly free combination Then, q e = q µ = 0 and q u = q d = q s = q c = q b = q τ /9. We are interested in Z ′ heavier than 4m µ coupled to aā such that its dominant decay mode is Z ′ → aā, leading to the µμµμ as in Eq. (2). Since in this option Z ′ does not couple to e and µ, the bounds from KLOE [20] and BABAR [17,21] will be relaxed. As Z ′ is heavier than 400 MeV and decays fast, bounds from the early universe or supernova cooling are irrelevant. In this case, q u can be even as large as O(0.1), leading to a large N production rate.
As mentioned above, ν µ produced in the pion decay will be in fact a linear combination of where U µ4 is the mixing between ν µ and N ′ (not with the heavier sterile particle, N ) so g νN can be as large as 10 −2 . On the other hand, the decay of heavier mesons such as K + or D + can produce both |ν 4 and |ν µ . The flux from the heavy meson decay can be decomposed as a flux of ν 4 proportional to |U µ4 | 2 plus a flux ofν µ proportional to (1 − |U µ4 | 2 ). While the scattering rate of ν 4 + quark → N + quark is proportional to g ′2 (1 − |U µ4 | 2 ), the scattering rate ofν µ + quark → N + quark should be proportional to g ′2 |U µ4 | 2 . Let us denote the ν µ flux predicted (in the absence of new physics) from the pion decay and the heavy meson decay respectively by F ν π and F ν H . The flux of N from the interaction of ν µ within the toy model described in the appendix should then be given by g 2 νN (F ν π + 2F ν H ). At FASERν for the ν µ energy lower than 1 TeV, we expect F ν π ≫ F ν H [1] so we will neglect the contribution from F ν H in our analysis. For ν µ energy higher than 1 TeV, the flux from the Kaon decay dominates, F ν π ≪ F ν H [1] so we will neglect F ν π , taking into account the extra factor of 2 in front of F ν H .

IV. DECAY OF NEW PARTICLES AND THEIR SIGNATURES AT FASERν
In this section, we first discuss the decay products of N assuming Z ′ is kinetically mixed with the photon. In the end, we shall comment on the case that Z ′ is the gauge boson of the B − 3L τ symmetry.
The m Z ′ > m N option: In the case m Z ′ > m N , the decay of N will be three-body with decay rate For m N > m N ′ + 2m µ , m N ′ + 2m π , the decay modes N → N ′ µ − µ + , N ′ π − π + also open up. In the absence of a faster decay mode for N (such as N → N ′ aā as discussed in the appendix), an N particle with an energy of E N can decay after traveling a distance of Considering that 0.2 m is smaller than the length of the detector, the majority of the decays take place within the detector in this parameter range. Since N is highly boosted in the beam direction, its track will be aligned with the beam direction so the signature will be a shower similar to that expected in the neutral current interactions of neutrinos and a pair of e − e + from a vertex separated from the shower vertex by ∼ 0.1 m. Since N is highly boosted in the forward direction, the line connecting the two vertices should make an angle of O(m N /E N ) or smaller with the beamline direction. The probability for one out of N N C ∼ 5000 background (SM) neutral current vertices being found behind the lepton vertex within a cone with an opening angle of m N /E N is less than (π/3)(m N /E N ) 2 (1.3/0.25) 2 N N C which will be less than one. That is despite the separation, we can identify which of the observed NC vertices is associated with a certain observed dilepton. If the kinematics allows, along the e − e + pair signal there will be also signals of µ − µ + or π + π − emitted close to the direction of the beam. However, the e − e + pair can come from the pair production by photons from the π 0 decay in the shower of neutral current interactions. Taking the crosssection of the pair production of ∼ 10 b/atoms [22], we find that the photons travel ∼ 1 cm before pair production. Moreover, while the invariant mass of e − e + from the photon will be of order of 2m e , that from N → N ′ e − e + will be of ∼ m N − m N ′ ≫ 2m e . Thus, by putting cuts on the distance between the shower vertex and the e − e + vertex and the invariant mass of e − e + , the background can be substantially reduced. Considering that the total neutral current events for FASERν are only O(5000), applying these cuts the number of background should become negligible.
Similarly, but with a rate suppressed by m 2 e /m 2 µ , the photons can produce muon pair [8]. The invariant mass of the muon antimuon pair from the photon will be close to 2m µ . As long as m N − m N ′ ≫ 2m µ , again by the measurement of the invariant mass of the charged lepton pair, the background can be vetoed. Another source of background for the µ − µ + pair is ν µ + nucleon → µ + c + X and the subsequent decay of c into µ. For contained events, this background can be reduced by the event topology (i.e., reconstructing the D meson track). Similar consideration applies for the background from ν e + nucleon → e + c + X and the subsequent decay c → e + X to the e − e + signal.
By measuring the distance between the neutral current vertex and the N decay vertex, one can derive some information on the N lifetime. However, since one of the final particles (N ′ ) will be missing, the derivation of the energy of N and hence the boost factor will not be possible. The derivation of the lifetime will also suffer from low statistics.
As discussed in sect. III and in the appendix, it is possible to obtain four muon signal through , both N and Z ′ should be heavier than 2m a (and hence heavier than 4m µ ) and moreover, the coupling of a to Z ′ should be much larger than eǫ which means N will decay promptly but the decay length of a can be within the range 1 mm−few meters. Again, the a andā particles will be highly boosted in the forward direction so the muon pairs will be emitted along the beam. If the event is fully contained in the FASERν detector by measuring the invariant mass of µ − µ + pair, the mass of a particle can in principle be derived. Moreover, by reconstructing the momenta of µ andμ, the directions of a andā and therefore the N decay vertex can be deduced. The distance between the N decay vertex and the a andā decay vertices is a measure of the a andā lifetime. However, for q ′ q = eǫq q smaller than the BABAR bound, the statistics of the contained event will be too low to perform such analysis. As discussed in the previous section, the four muon events originated in the rock before the detector can increase the statistics. We will discuss this possibility further in sect V.
where f can be e or µ and the factor of 3 in the denominator comes from averaging over three polarizations of the initial vector boson. Z ′ can also decay into π + π − but we shall focus only on the leptonic decay modes. The decay length of Z ′ is of order of 0.1 cm MeV, the decay length will be shorter than 1 mm, so FASERν cannot disentangle the track of Z ′ . In this case, Z ′ can decay into µ + µ − so if the ν µ scattering takes place inside the detector, the signal will be a nuclear shower and a µ − µ + pair from a single vertex which is background free. By measuring the energy momentum of µ + µ − , m Z ′ can be reconstructed. Another signal will be a nuclear shower plus e − e + pair. Again the invariant mass of the e − e + pair gives m Z ′ . Since Z ′ will be highly boosted, the angle between the final fermion pair will be small and of order of m Z ′ /E Z ′ ∼ 10 −3 but the angular resolution of FASERν is better than 0.1 milliradian [1] so this angle can be reconstructed. We can write where θ is the angle between f andf . Since there will be an uncertainty of 30% in the reconstruction of p f and pf (i.e., δp f /p f ∼ 30% [1]), the uncertainty of m Z ′ derived from a single pair of ff will be 40%. Of course, if the statistics is high, the uncertainty in the derivation of m Z ′ will be reduced. If N pair pairs are registered, the uncertainties will be 40%/ N pair . As discussed in the appendix and above, by coupling the Z ′ particle to a pair of a scalars lighter than m Z ′ /2, we can have a dominant four-muon signal which is background free. This of course requires Z ′ to be heavier than 4m µ .
Notice that in both cases m Z ′ > m N and m N > m Z ′ , the main decay mode of N produces N ′ which is a metastable particle. As discussed in sect. III, the lifetime of N ′ will be long enough to exit the detector so N ′ will appear as missing energy-momentum.
As mentioned in sect. II, the neutrinos can interact in the soil and concrete before reaching the detector. In this case, the e − e + or π − π + pairs (as well as the neutral current showers) will be absorbed but the muon pairs can reach the detector. Since the angular deviation of muon and antimuon during the propagation will be much larger than the angle that they make with each other at the production (m Z ′ /E Z ′ ∼ 10 −3 ), these through-going muons cannot be used to derive the mass of Z ′ ; however, they will be indicative of new physics and taking them into account a stronger bound on the coupling can be obtained. We shall study this in more detail.
Let us now briefly discuss the case that Z ′ is the gauge boson of B − 3L τ . If the kinematics allows Z ′ can decay into π + π − or τ + τ − but as discussed before we are interested in the case that the decay mode N → N ′ Z ′ ( * ) → N ′ aā dominates and we obtain a four muon signature. The rest of the discussion is similar to above, except that here the BABAR bounds do not apply so the statistics of the contained vertex events can be large enough to derive information on the a andā lifetime through the method that we discussed before. By measurement of the invariant mass of the contained vertex muon pairs, the mass of a can also be derived with a precision of 40%/ N pair where N pair is the number of contained vertex muon pair. For m Z ′ < m N , the invariant mass of µμµμ gives the mass of m Z ′ with a precision of (δE µ /E µ ) 4/N pair = 60%/ N pair . For m Z ′ > m N , the invariant mass of µμµμ will have a continuous spectrum and will not be peaked at m Z ′ but m a can still be derived by extracting the invariant mass of the µμ pairs.

V. PRODUCTION OF N VIA NEUTRINO FLUX AND THE SIGNAL AT FASERν
In this section, we first discuss the N production rate assuming Z ′ is kinetically mixed with the photon. The same discussion holds valid for the case that Z ′ is the gauge boson of the B − 3L τ symmetry, replacing eǫq q with the common gauge coupling of the quarks. The cross section of the scattering of ν µ ≃ν µ = 3 i=1 U µi ν i off quarks can be written as where g νN = g ′ U µ4 and s is the Mandelstam variable for the system of ν µ and parton carrying a fraction x of the proton momentum, s = 2xm p E νµ ∼ (10 GeV) 2 (x/0.1)(E νµ /500 GeV) and θ is the scattering angle in the center of mass frame of the ν µ -parton system. q q is the electric charge of the parton; i.e., q u = 2/3 and . To obtain the ν µ cross section off the nucleus, we should convolute with the Parton Distribution Functions (PDFs) where F q and Fq are parton distribution functions and t = (m 2 N − s)(1 − cos θ). Notice that as θ → 0, t goes to zero. This can be understood as we have neglected the masses of the partons and the t variable associated with the scattering of a massless particle to another in the forward direction vanishes. In the limit t → 0, dσ/d cos θ in Eq. (5) will be dramatically enhanced. In fact, the total cross-section is , in which s = 2xm p E νµ . The 1/m 2 Z ′ behavior of σ(ν µ + q → N + q) shows that for the majority of the scatterings, |t| ∼ m 2 Z ′ or smaller. On the other hand, the dependence of PDFs on t is mild. We can therefore simplify the integration in Eq. (6) by setting F q (x, t) = F q (x, −m 2 Z ′ ) and Fq(x, t) = Fq(x, −m 2 Z ′ ) and write Notice that we have used the fact that the cross sections of the scattering off quark and antiquarks are equal. Similarly, the scattering cross sections of ν µ andν µ off quarks are equal. The number of ν µ andν µ converting to N inside FASERν can therefore be estimated as M det is the total mass of FASERν. F νµ (E µ ) and Fν µ (E µ ) are the fluxes of ν µ andν µ at the detector (integrated over time). F νµ (E µ ) and Fν µ (E µ ) can be read from Fig. 4 of [1]. For E ν <TeV, the majority of the neutrino flux comes from the pion decay. Since pion is lighter than m µ + m N ′ , the flux will be purely composed ofν µ ≃ i U µi ν i . For E ν >TeV, the Kaon decay dominates so the flux at the detector can be decomposed as theν µ part which takes a fraction of (1 − |U µ4 | 2 ) of the flux plus the N ′ part which constitutes a fraction of |U µ4 | 2 of it. Since the cross section of N ′ is |U µ4 | −2 times that ofν µ , for E ν >TeV, the contribution from N ′ to N HESE should be equal to that fromν µ . The Θ-function in Eq. (8) accounts for the contribution from N ′ .
Let us now study how much the discovery reach of FASERν can be improved by including ν µ + nucleon → N + X taking place in the rock before the detector. As discussed before for the through-going events, we cannot discriminate the signal and background dimuon events. We therefore assume m N > m N ′ + 4m µ and suppose the signal is the background-free µμµμ as described in the appendix. As mentioned before, there are 10 meters of concrete and 90 meters of rock before the detector. As demonstrated in the right panel of Fig. 5 of [1], the ν µ beam reaching the FASERν detector is very well collimated along the proton beamline. For simplicity, we shall assume that all ν µ flux is directed towards the center of the detector. At the scattering of ν µ off a parton, the angle in the lab frame θ lab (the angle between the produced N and the initial ν µ ) will be smaller than m Z ′ /E ν < 10 −3 . The angular spread of the muon particle from the direction of N can be estimated as θ ′ lab ∼ (m N − 3m µ )/(2E N ) < few × 10 −3 . As mentioned before the angular dispersion due to propagation of the produced muons is of order of 0.5 • = 0.008 radians [7] which is larger than θ lab and θ ′ lab . As we discussed before, in order to produce two muon pair signals (µμµμ), the Z ′ particle has to have a large coupling to an intermediate scalar, a. Thus, the whole N → N ′ Z ′ ( * ) → N ′ aā process will take place promptly. However, the a particles can travel a sizable distance before producing µμ. If the production of the µμ pair takes place at a distance smaller than (0.25 m/2)/0.008 = 14 m, practically both µ − and µ + from the a decay will reach the detector. If the production takes place farther, we should multiply the rate with the probability that both µ − and µ + reach the detector. If the production takes place at a distance z > 14 m from the detector, the probability of each of µ − or µ + arriving at the detector will be so the probability for both muons reaching the detector will be given by p 2 µ . Remember that if only one of them reaches the detector, it cannot be distinguished from the background muons originating from the IP. In fact, both pairs from the a andā decay have to pass through the detector in order to distinguish the signal event from the background. Let us suppose a andā respectively decay at distances of z 1 and z 2 from the detector. The probability that all the four muon and antimuons arrive at the detector is given by [p(z 1 )] 2 × [p(z 2 )] 2 .
As long as the distance traveled by a andā before decay is shorter than 10 m, the number of muon pairs reaching the detector originated from the en-route rock and concrete can be estimated as where we have assumed that practically all N decays lead to a µμµμ. ρ soil and ρ con are respectively the soil and concrete densities. To carry out our computation, we take ρ soil = ρ con = 2.5 gr/cm 3 . As discussed in sect II, if FASERν is equipped with the scintillator plate in its front and records the timing of the arrival of the incoming µμµμ event, the background will be negligible. As a result, even the detection of a single µμµμ event can be regarded as an indication for new physics. In fact, even if only three out of four µμµμ reaches the detector, the signal will be background free. The sum of the numbers of the µμµ,μµμ and µμµμ events is given by Eq (10), replacing [p µ (z)] 4 with [p µ (z)] 3 . In our analysis, we will however focus only on four muon events.

VI. RESULTS
Let us now study the number of events at FASERν and then discuss what bounds can be derived from this experiment. For the time integrated neutrino fluxes, we have used Fig 4 of [1] which corresponds to the run III of the LHC during 2021-2023. For Parton Distribution Function (PDF), we have taken the CT10 model [23] from LHAPDF-6.3.0 software [24].   1 shows the number of produced N versus m Z ′ for 20 MeV < m Z ′ < 10 GeV, taking g νN = 10 −2 and ǫ = 7 × 10 −4 . To draw this plot, we have set m N = 5 GeV so, as long as m Z ′ > 4m µ the kinematics allows the production of four muons from N decay. As seen from these figures, including the through-going events increases the statistics by a factor of 4. For m Z ′ < 1 GeV (0.5 GeV), N HESE + N through−going (N HESE ) will exceed 1 so hints for new physics will be revealed in Run III of the LHC. Notice that m N = 5 GeV is well beyond the reach of MINERνA and CHARM II experiments [25]. Fig. 2 shows the number of produced N versus m N for various values of m Z ′ , again taking g νN = 10 −2 and ǫ = 7 × 10 −4 . As seen from this figure, if m Z ′ ∼ 0.1 GeV, m N up to a few 10 GeV can be probed by FASERν during run III of the LHC so FASERν will give unprecedented bounds. For m Z ′ >few GeV, the number of events will be less than 1 so FASERν during the run III cannot test Z ′ heavier than ∼ 2 GeV but during the high luminosity run of the LHC and with an upgrade of FASERν, the parameter range with heavier values of Z ′ can also be probed. For m N < 4m µ , the kinematics does not allow the background free µμµμ signal so the through-going events cannot help to extend the effective volume of FASERν. For m N < 2m µ , the µμ signal will also be absent but e − e + can be produced and detected inside FASERν.
The blue and red dots in Fig. 2 correspond to the best fit of the MiniBooNE solutions, respectively, in Ref. [5] with ǫg νN = 2.7 × 10 −6 , m N = 150 MeV, m Z ′ = 1.2 GeV and in Ref. [6] with ǫg νN = 2 × 10 −8 , m N = 130 MeV, m Z ′ = 30 MeV. Notice that at MiniBooNE, the neutrino energy in the beam is much lower. As a result of light Z ′ [6], the scattering off nucleons in the nucleus at MiniBooNE is coherent, leading to a huge enhancement of the cross-section by the square of the atomic number for a given coupling. At FASERν, the energy-momentum exchange will be too high to maintain coherence so we shall not have such enhancement. As seen from the figure, the number of events during run III of the LHC at FASERν for the values of parameters providing a solution to MiniBooNE will be less than 0.1. With an upgrade of FASERν during the high luminosity run of the LHC, the statistics can reach 1000 times larger [1]. These solutions can then be tested.
As shown in the figure, as m N varies between 0.1 GeV to 6 GeV, the number of events only slightly changes. In general, for s ≫ m 2 N , we expect only a weak dependence on m N but for s ∼ m 2 N , the dependence should be strong so the weak dependence of the number of events on m N means the main contribution to σ tot νN comes from relatively large values of x for which s ≫ m 2 N . As seen from Fig. 1, the dependence of the number of events on m Z ′ is however strong. In the limit s ≫ m 2 Z ′ , from Eq. (7), we find σ ∝ 1/m 2 Z ′ which receives the dominant contribution from small scattering angles; i.e., θ < ∼ m Z ′ / √ s in Eq. (5) where θ is the scattering angle in the center of mass frame of the parton neutrino system. The corresponding scattering angle in the lab frame is θ lab < θ √ s/E ν = m Z ′ /E ν = 10 −3 (m Z ′ /1 GeV). The dashed (solid) red curves show the reach of FASERν during run III of the LHC, including (without) the through-going signal. That is to draw the red solid (dashed) curve, we have set N HESE = 1 (N HESE + N through−going = 1). The blue curves marked with FASER2ν show the improvement on the bounds if the data increases by a factor of 1000. Such an increase is feasible during the high luminosity phase of the LHC as the integrated luminosity will increase twenty times in this phase of the LHC and the mass of the detector is under discussion to be increased by a factor of ten to thousand times [1]. The results are derived under the assumption of zero background.
The dashed horizontal black lines in Fig. 3 show the present combined constraint on ǫ from BABAR [17] and on g νN from theoretical consideration (see the appendix): ǫg νN < 7 × 10 −6 . As seen from the figures, for m N < few 10 GeV, FASERν can probe the values of the coupling well below this combined bound. As discussed before, for m N <GeV, the scenario could also be tested by MINERνA and CHARM II (as shown by the cyan and magenta solid curves taken from [25]) but the range 2 GeV < m N will be explored by FASERν for the first time. The green region shows the solution to the MiniBooNE anomaly [25]. Upgrade of FASERν with 1000 times more data can probe the entire this region. Notice in the model that the coupling of the Z ′ with the standard model fermions is through a kinetic mixing with the photon, q ′ u = −2q ′ d = −2q ′ s = (2eǫ/3) but in the gauged B − 3L τ symmetry, q ′ u = q ′ d = q ′ s . However, for m N , m Z ′ > 4m µ when the µμµμ signal mode is open, the bounds on g νN q ′ q from FASERν for both models will be very similar. In the latter case, since the present combined bound on g νN q ′ q is relatively weak and of order of 10 −3 , FASERν during run III of the LHC can improve the bound by two orders of magnitude. IF FASERν finds a large number of µμµμ signal exceeding few 100 according to Figs. 1 and 2, g νN q ′ q should be much larger than 7 × 10 −6 so the option of Z ′ mixed with the photon will be ruled out, providing a hint in favor of the option of Z ′ as the gauge boson of the local B − 3L τ symmetry.  [25]. The green area shows the solution to the MiniBOONE anomaly presented in [6]. The dashed black lines show 7 × 10 −6 which is a combination of the bound ǫ < 7 × 10 −4 from BABAR [17] and setting gνN = 10 −2 .

VII. SUMMARY AND CONCLUSIONS
We have studied the discovery potential of FASERν for beyond standard model interaction of neutrinos with nuclei that leads to a multilepton signature. In this model, the interaction of the neutrino flux creates a heavier fermion, N whose decay produces the multi-lepton signal. Similarly to the models proposed in [5,6,9], the ν + nucleon → N + X scattering takes place through the exchange of a new U (1) gauge boson, Z ′ with mass smaller than ∼ 1 GeV and kinetically mixed with the photon. The decay of the produced N can then produce lepton pair via either on-shell (for m Z ′ < m N ) or off-shell (for m Z ′ > m N ) Z ′ .
We have shown how to build a consistent model in which the coupling between ν, N and Z ′ is relatively large (of order of 10 −2 ), respecting all the present bounds. With such a coupling, we have found that the FASERν detector can record between 1000 to 1 events as m Z ′ varies between 10 MeV to 1 GeV. As long as m N < ∼ 10 GeV, the dependence of the number of events on m N is only mild because at FASERν, the center of mass energy in the scattering process is much larger than the N mass. The part of the parameter space of the model with m N < 1 GeV has already been probed by CHARM II, MINERνA and MiniBooNE experiments but heavier N can be probed by FASERν for the first time. We argue that the dilepton signal (i.e., the e − e + or µ − µ + signals coming from a single vertex) produced by the neutrino flux interaction inside FASERν can be discriminated against the background from the pair production by photons and from charmed-induced events (i.e., ν µ (ν e ) → µ(e) + c + X and the subsequent c decay into µ (e)), respectively, by reconstructing the invariant mass of the lepton pair and by determining the event topology. After applying the cuts, the signal will be practically background free. The bound that FASERν can set on the coupling is shown in Fig 3. The neutrino beam before reaching the detector has to pass through 100 meters of rock and concrete. For the first time, we discussed the possibility of enlarging the effective mass of the detector to search for new physics by including the events originated in the rock. All the charged particles except for muons will be absorbed in the rock so we focused on the multimuon signal. Unlike the case of dimuon events originating inside the detector, the dimuon events from the rock will suffer from a large background. Equipping the front side of the detector with a scintillator plate capable of recording the timing of the event can significantly reduce the background; however, the background for the through-going dimuon events will be still large. Unlike the case of the events starting inside the FASERν detector, we cannot employ cuts on the invariant mass of the µ − µ + pair or on the event topology to reduce the background for the dimuon signals. However, the multi-muon events with more than two muons will become background free. Fig 3 shows how the bounds can be improved by including the through-going events in the search for the four-muon signal.
We showed how by adding a new scalar to the model a µμµμ signal can be achieved. In this variation of the model, the coupling of Z ′ to the scalar dominates over the coupling to µ + µ − so the N decay dominantly produces a pair of these scalars rather than the µμ pair. The produced a andā eventually decay into µμ comprising the µμµμ signal.
Throughout the paper, we have focused on the scattering of the ν µ flux. In principle, ν e and ν τ can also have large couplings to N and Z ′ . Similar arguments can be repeated for the ν e and ν τ scattering, too. However, since the flux of ν e (ν τ ) is smaller than that of ν µ by a factor of 15 (1000), the bounds that can be derived on the relevant coupling will be weaker by a factor of about 4 (30) compared to the bound on the ν µ coupling to Z ′ .
We have also discussed the option that Z ′ is the gauge boson of a B − 3L τ local symmetry with negligible coupling to µ and the electron. In this case, the bounds from BABAR on the Z ′ coupling to the SM fermion is relaxed so the ν µ scattering cross-section off nuclei can be large enough to produce more than 20000 N even when both m Z ′ and m N are heavier than 4m µ . If the observed number of µμµμ signal events by FASERν during the run III of the LHC exceeds a hundred, this option should be taken more seriously as for the case of Z ′ mixed with the photon, the number of µμµμ events cannot be such large. We have shown that if the statistics are enough, the parameters of the model such as the lifetime of a and its mass as well as the mass of Z ′ can be extracted from the data. In the case of the null signal at FASERν, the bound on g νN q ′ q can be improved by two orders of magnitude. The idea of using the through-going muons to enlarge the effective mass of the detector is quite general and can be applied in various contexts beyond the four-muon signal produced by neutrino scattering. For example, if either of a andā decays into e − µ + (instead of into µ − µ + ), we can have three-muon through-going signal which is again background-free. The multi-muon through-going signal may originate from the decay of a new particle produced at IP rather than by neutrino interaction in the rock. Finally, studying the through-going dimuon events will increase the ν µ + nucleon → µ + c + X data sample. Of course, eliminating the background from the pile-up of the muons from IP requires the timing of the arrival of the muons comprising a multi-muon signal which in turn motivates installing a scintillator plate in front of the detector.

APPENDIX
In this section, we introduce a model for the interaction of form shown in Eq. (1). The idea is based on the model which was introduced in Ref. [26]. In this model, there are two left-handed sterile neutrinos N L and N ′ L with an off-diagonal coupling to the gauge boson of the new U (1) gauge symmetry, Z ′ , as follows This form of the interaction can be easily obtained from the gauge symmetry if we simply assign opposite charges to If N ′ L mixes with ν µ , an interaction of form (1) can be achieved with g νN = g ′ U µ4 . Notice that we do not want N L to mix with active neutrinos as there are strong bounds on such mixing for N with a mass of few 100 MeV-few GeV. As we discussed in sect III as long as 30 MeV < m N ′ < 70 MeV, the bounds on the mixing with ν µ is much more relaxed. The difference in mixing of N L and N ′ L requires breaking of the U (1) gauge symmetry. Let us introduce a new scalar, φ which is neutral under the standard model gauge group but, under the new U (1) gauge symmetry, has a charge equal to that of ψ 2 . Moreover, let us impose a Z 2 symmetry under which As a result, the combinations φψ 1 − φ * ψ 2 and φψ 1 + φ * ψ 2 are respectively even and odd under the Z 2 symmetry. To give a Dirac mass to N L and N ′ L , we add N R and N ′ R which both are singlets of the gauge groups but have opposite Z 2 parities: N ′ R (N R ) is Z 2 even (odd). We can then write the following Yukawa couplings that preserve both the gauge symmetry and the Z 2 symmetry: Without loss of generality, we can invoke the global U (1) symmetry to rephase φ and make its vacuum expectation value, φ = v φ / √ 2, real. These terms then lead to Dirac mass terms as follows where Since N ′ R is Z 2 even, we can write the following Yukawa coupling: where H is the standard model Higgs, L α is the left-handed lepton doublet of flavor α. Since we are mostly interested in ν µ , we may identify α with µ. The Z 2 symmetry forbids writing a similar Yukawa coupling for N R so only N ′ mixes with ν α . The induced mass term can be written as This mass matrix, which has a zero mass eigenstate (which is mainly composed of the active neutrino, ν α ), can be diagonalized by Thus, the mixing of ν α and N ′ will be given by U α4 = Y Rα H /m N ′ . Notice that, unlike the minimal model when N ′ L directly couples to ν α , in our model where the coupling is through N ′ R , the mixing does not lead to a contribution of m N ′ U 2 α4 to the ν α mass so the mixing can be relatively large (U α4 ∼ 0.1). Notice that the gauge and Z 2 symmetries allow a mass term of form This mass term should appear as the (2, 2) element of the mass matrix in Eq. (16) and would break the lepton number, inducing a Majorana mass for ν α proportional to µ like in the inverse seesaw mechanism [27].
The N ′ particles with a mass of 50 MeV and a mixing of |U µ4 | ∼ 10 −2 can be produced in the supernova core. Since they will reach equilibrium with the matter, the bounds from supernova cooling do not rule out this model. Similarly, the N ′ abundance in the early universe will be sufficiently reduced by scattering N ′ + f → ν µ + f when the temperature drops its mass. As a result, this part of the parameter space is not restricted by cosmological data [15].
Remember that in order to have large signal sample at FASERν, we prefer the parameter range m N ′ ∼ 50 MeV, m N ∼ 1 GeV and m Z ′ /g ′ ∼ 50 MeV. To obtain m N = Y v φ / √ 2 ∼ 1 GeV, we need v φ > ∼GeV. m N ′ = Y ′ v φ / √ 2 ∼ 50 MeV can then easily be obtained with Y ′ /Y ∼ 0.01. However, m 2 Z ′ also receives a contribution given by g ′2 v 2 φ /2. With m Z ′ /g ′ ∼ 500 MeV, we can still have a handful of signal events at FASERν (see Fig. 1). To have m Z ′ /g ′ ∼ 50 MeV (and therefore to obtain a large number of events), there should be a cancellation. Of course, vacuum expectation values from additional scalars charged under new U (1) will only increase m 2 Z ′ as the contribution from each will be positive. However, if we invoke the Stuckelberg mechanism, we can obtain a negative contribution to m 2 Z ′ which cancels out the positive contribution from v 2 φ . The smaller m Z ′ , the higher degree of fine tuning is required. The range m Z ′ /g ′ ∼ 500 MeV seems to be more natural. Moreover, lighter Z ′ can only decay to the e − e + pair not being able to produce multi-muon events.
If φ is lighter than m Z ′ /2, Z ′ can decay into φφ faster than into the lepton pair because g ′ ≫ eq ′ f . φ can promptly decay intoN ′ R N ′ L via the relatively large Y ′ coupling. In principle, after the electroweak and the new U (1) symmetry breaking through a coupling of form λ φH |φ| 2 |H| 2 , φ can mix with the Higgs and therefore obtain a coupling of form λ φµ φμµ. The coupling will be however suppressed by λ φµ ∼ λ φH v φ m µ /m 2 H = 10 −5 λ φH ≪ Y ′ so the dominant decay mode of φ will be invisible φ →N ′ R N ′ L . To avoid missing the signal, we can take φ to be heavier than m Z ′ /2 so that the dominant decay mode of Z ′ would be decay to the lepton pair. If we want the Z ′ decay to produce µμµμ instead of just one pair of µμ, we may introduce another scalar, a, singlet under gauge symmetries, lighter than m Z ′ /2 and mixed both with the Higgs and with φ through A φa a|φ| 2 + A Ha a|H| 2 + H.c.
The mixings of a with φ and H can approximately be written as α ∼ A φa v φ /m 2 φ and β ∼ A Ha v/m 2 H , respectively. As long as g ′ sin α > q ′ f , the dominant decay mode of Z ′ will be decay into an a pair. The a particle will have a Yukawa coupling of form ( √ 2m µ sin β/v)aμµ. As long as sin β > ∼ 10 −3 , the decay length of a with an energy of few 100 GeV will be smaller than ∼ 10 m. Moreover, for sin β < ∼ 10 −3 (|A Ha | < ∼ 10 −3 m H ), the A Ha coupling does not induce any instabilities for the Higgs potential. If we want to obtain larger coupling between a and the µμ pair, we can invoke a heavy inert Higgs doublet, H I with large coupling to the muons and a trilinear coupling of form aH † I H.