Discovering Dark Matter at the LHC through Its Nuclear Scattering in Far-Forward Emulsion and Liquid Argon Detectors

The LHC may produce light, weakly-interacting particles that decay to dark matter, creating an intense and highly collimated beam of dark matter particles in the far-forward direction. We investigate the prospects for detecting this dark matter in two far-forward detectors proposed for a future Forward Physics Facility: FASER$\nu$2, a 10-tonne emulsion detector, and FLArE, a 10- to 100-tonne LArTPC. We focus here on nuclear scattering, including elastic scattering, resonant pion production, and deep inelastic scattering, and devise cuts that efficiently remove the neutrino-induced background. In the invisibly-decaying dark photon scenario, DM-nuclear scattering probes new parameter space for dark matter masses 5 MeV $\lesssim m_{\chi} \lesssim$ 500 MeV. When combined with the DM-electron scattering studied previously, FASER$\nu$2 and FLArE will be able to discover dark matter in a large swath of the cosmologically-favored parameter space with MeV $\lesssim m_{\chi} \lesssim $ GeV.

A primary goal of high-energy colliders is to produce dark matter (DM) particles. If DM is heavy with a mass near the weak scale, its signature is missing transverse energy, which has been studied in detail for decades. If DM is light, however, such searches are typically ineffective (as are conventional direct detection searches), and alternative search strategies, experiments, and facilities are needed.
In this study, we consider extremely simple models of light DM in which the Standard Model (SM) is supplemented by a dark photon [1] that decays to pairs of DM particles through A → χχ. For dark photons with typical loop-suppressed couplings ε ∼ 10 −4 − 10 −3 and m A , m χ ∼ MeV − GeV, the DM annihilates through χχ → A ( * ) → ff in the early universe, yielding the correct thermal relic density. This model is representative of a broad class of hidden sector theories in which the correct amount of DM is produced through thermal freeze-out within the standard cosmology [2][3][4][5][6][7], just as in the case of weak-scale DM. In this scenario, however, the DM is light. As a result, at colliders, the dark photons and DM are dominantly produced along the beampipe in the far-forward region, escape through holes in collider detectors, and evade all conventional collider searches.
To remove such "blind spots" from the Large Hadron Collider (LHC) physics program, a number of experiments are currently planned for the far-forward region. FASER [8][9][10][11] has been completely constructed, and FASERν [12][13][14] and SND@LHC [15] are also being prepared to take data when Run 3 of the LHC begins in 2022. For the High Luminosity-LHC (HL-LHC) era, a Forward Physics Facility (FPF) is under study [16][17][18]. The FPF would house a suite of far-forward experiments, including possibly FASER2 [19], targeting new long-lived particles that decay visibly in the detector; FORMOSA [20], a milli-charged particle detector; FASERν2 [21,22], a 10-tonne emulsion detector; SND2, a successor to SND@LHC; and FLArE [23], a proposed liquid argon time projection chamber (LArTPC) with an active volume of 10 tonnes (FLArE-10) to 100 tonnes (FLArE-100). FASERν2, SND2, and FLArE will detect millions of TeV-energy neutrinos, providing a wealth of SM measurements, but they also have the potential to search for light DM and other new particles.
Here we evaluate the prospects for discovering light DM at FASERν2 and FLArE through DMnuclear scattering in the HL-LHC era. This work complements Ref. [23], which focused on the prospects for observing elastic DM-electron interactions in these detectors; Refs. [24,25], which explored the potential of FASER to probe inelastic DM; Ref. [26], which studied the scatterings of unstable, but very long-lived, heavy neutral leptons at FASERν2; and Ref. [27], which investigated leptophobic DM scattering at SND@LHC. 1 We assume these experiments are located in a new cavern that is under study for the FPF, which would place the fronts of these detectors approximately 620 m from the ATLAS interaction point (IP), and we consider 14 TeV pp collisions and the expected HL-LHC integrated luminosity of 3 ab −1 . Alternative locations for the FPF that are ∼ 150 m closer or farther from the IP do not change the prospects much, provided, of course, that they are large enough to house the detectors we consider.
We begin by defining the light DM models in Sec. II and specifying the detectors in Sec. III. We then consider the dominant processes contributing to DM-nuclear scattering, including elastic scattering (χp → χp), resonant pion production (χN → χN π), and deep inelastic scattering (DIS) (χN → χX) in Secs. IV, V, and VI, respectively. For each of these signals, we devise simple kinematic cuts to differentiate the DM signal from the neutrino-induced SM background.
In Sec. VII, we then combine all of these DM-nuclear probes with the DM-electron signals investigated in Ref. [23]. We find that DM-nuclear scattering and DM-electron scattering are quite complementary, with nuclear scattering more powerful for relatively high masses m χ 100 MeV and electron scattering more sensitive for low masses m χ 10 MeV. By combining DM-nuclear and DM-electron scattering, FASERν2 and FLArE can cover the cosmologically-favored parameter space, where the χ thermal relic density is at or below Ω DM , for a wide range of DM masses between MeV m χ GeV. In Sec. VII, we also compare the sensitivity of FASERν2 and FLArE to non-LHC experiments that have discovery potential for invisibly decaying dark photons and light DM [42,43]. Our conclusions are presented in Sec. VIII.

II. INVISIBLY-DECAYING DARK PHOTON MODELS
In this section, we describe two popular benchmark models in which light DM interacts with the SM through an invisibly decaying dark photon mediator. Given its coupling to electrically charged particles and quarks, in particular, the dark photon efficiently mediates scattering between DM and nuclei, making these models an interesting test case for our study.
The dark photon, A , is a massive gauge boson that arises when the SM is supplemented with a new broken U(1) D symmetry. For light GeV-scale dark photons, the dark photon Lagrangian is where F µν is the dark photon's field strength, m A is the dark photon mass, ε is the kinetic mixing parameter, J µ EM and J µ D are the SM electromagnetic and U(1) D currents, respectively, and g D ≡ √ 4πα D is the U(1) D gauge coupling. For the DM candidates, χ, we will examine two popular examples: Majorana fermion DM and complex scalar DM. The corresponding Lagrangians are where m χ is the DM mass. The U (1) D currents associated with these models are These two DM models have many similarities, but also some key differences. We discuss them in turn, beginning with the Majorana fermion case. As noted in Sec. I, an attractive feature of these light DM models is the fact that the observed DM relic density can be easily obtained through thermal freeze-out. For m A > 2m χ , Majorana fermion DM annihilates in the early universe through χχ → A ( * ) → ff with cross section where we have assumed m A m χ and y ≡ ε 2 α D (m χ /m A ) 4 [5]. As evident from Eq. (4), the annihilation is P -wave, and so bounds from cosmic microwave background (CMB) temperature anisotropies on late-time DM annihilation are not very constraining in these models [44,45]. In addition, the scattering of Majorana fermion DM in direct detection experiments is also velocitysuppressed at the non-relativistic energies relevant for these searches, and so direct detection null results also do not set strong limits.
For complex scalar DM, the annihilation cross section is, in fact, similar to that for Majorana fermion DM. Equation (4) still applies, and so the complex scalar DM model also evades CMB bounds. In contrast to the Majorana fermion case, however, the non-relativistic scattering of complex scalar DM in matter is not velocity-suppressed. Direct detection null results are therefore a significant constraint on this model. These bounds may be evaded, however, if a small mass splitting is introduced to make the DM scattering transition inelastic [46].
In this work, we will present our results in the (m χ , y) plane. As we will see, at the relativistic energies relevant for the LHC, the DM-nuclear interactions for Majorana fermion and complex scalar DM are very similar, and so the results we derive will be almost imperceptibly different in the (m χ , y) plane. We will therefore simply present the Majorana fermion DM results. At the same time, to understand the cosmological significance of these results, we will also present "thermal targets," the regions of parameter space where the thermal relic density is identical to the observed DM abundance. These differ slightly for the Majorana fermion and complex scalar DM models, and so we will present both, using the relic density predictions of Ref. [47].
To reduce the parameter space to two dimensions, we will present results for α D = 0.5 and m A = 3m χ throughout this work. These represent relatively conservative choices in terms of characterizing the experimental prospects for testing the thermal freeze-out hypothesis, at least in the regime m A m χ . Of course, if m A − 2m χ m A , the annihilation rate is resonantly enhanced, and the corresponding thermal targets occur at smaller couplings and can be much more challenging to probe at colliders [48][49][50].

A. Benchmark Detectors
The benchmark detectors we consider are identical to those studied in Ref. [23], except that they are now assumed to be housed in the "new cavern" FPF, placing them 620 m from the ATLAS IP. We review their most salient characteristics here; for more details, see Ref. [23].
FASERν2 [21] is envisioned to be a larger version of FASERν [12], currently being built for LHC Run 3. The FASERν2 benchmark detector we consider here is a 10-tonne rectangular tungstenemulsion detector with location and size given by where L is the distance from the IP to the front of the detector, and ∆ and S T are the longitudinal and transverse dimensions of the tungsten target. At the ATLAS IP during the HL-LHC, it is expected that the beam half-crossing angle will vary by as much as 250 µrad, moving the beam collision axis horizontally by as much as 15 cm at the detector location. Given the detector's transverse dimensions and the ∼ 20 cm spread of the DM signal and neutrino background [51], the crossing angle will have little effect on our results; for simplicity, we assume that the detector is always centered on the beam collision axis. We will assume that tracks down to momenta of 300 MeV can be detected and that the emulsion is exchanged periodically so that the track density remains manageable. This requires changing the detector every 30 fb −1 or so, or less if a sweeper magnet is available to bend away muons produced at the IP.
The main disadvantage of emulsion detectors for this DM search is the lack of timing, which makes it difficult to reject muon-induced backgrounds. To remedy this, it is necessary to augment the tungsten-emulsion detector with interleaved electronic tracker layers, which would provide event time information. This design follows the successful example of the OPERA experiment [52], and an analogous design is being implemented for SND@LHC [53]. We will, therefore, assume that muon-induced backgrounds can be rejected by vetoing events in coincidence with a high-energy muon track. It is important to note, however, that all of our FASERν2 sensitivities depend on this assumption, and if muon-induced backgrounds are difficult to reject in emulsion detectors, liquid argon technology may be preferable for dark matter detection.
For FLArE, we consider two sizes with physical dimensions where, as above, L is the distance from the IP to the front of the detector, ∆ and S T are the longitudinal and transverse dimensions of the active volume, and we assume that the detector is centered on the beam collision axis. Particle kinetic energy thresholds for LArTPC detectors typically lie in the 10-100 MeV range. For protons, we will consider two kinetic energy thresholds: a conservative value of 50 MeV, as is considered in the DUNE Conceptual Design Report [54], and a more optimistic choice of 20 MeV. Concerning the latter, we note that the ArgoNeuT experiment has already achieved thresholds for such short proton tracks down to 21 MeV [55,56]. For other particles, including shower-like objects (electrons, photons, neutral pions) and charged pions, we will assume a 30 MeV kinetic energy threshold, which is broadly consistent with Refs. [54][55][56]. In contrast to emulsion detectors, LArTPCs have good active event timing capabilities, particularly when equipped with a light collection system [57,58], and we will assume that vetoing events with a coincident muon is sufficient to remove all muon-induced backgrounds.

B. Expected Neutrino Fluxes
A crucial ingredient for the estimation of background rates is the flux of neutrinos passing through the different detectors. We use the dedicated forward physics event generator Sibyll 2.3c [59][60][61], as implemented in the CRMC simulation package [62], to simulate the primary collisions. We then use the fast neutrino flux simulation introduced in Ref. [51] to simulate the propagation of SM hadrons through the LHC beam pipe and magnets and their decays into neutrinos.
The results are presented in Fig. 1 for the HL-LHC with an integrated luminosity of 3 ab −1 , assuming no beam crossing angle. The upper panels show the numbers of neutrinos passing through the detectors. Unsurprisingly, detectors with a larger cross sectional area will have more neutrinos passing through them. The lower panels show the numbers of charged current (CC) and neutral current (NC) DIS neutrino interactions in the detectors, where we use the neutrino interaction cross sections from Ref. [13]. Note that the event rate is larger for FASERν2 than FLArE-10, despite the two detectors having the same mass. This is because the neutrino beam is strongly collimated around the beam collision axis, and so a narrow detector with more mass close to the beam axis, such as FASERν2, will observe a larger event rate. During the HL-LHC era, we expect about 3.9 × 10 4 electron neutrino, 2.2 × 10 5 muon neutrino, and 1.5 × 10 3 tau neutrino CC interactions in the FLArE-10 detector. In addition, we expect about 8.9 × 10 4 NC neutrino interactions. The average energy of these interacting neutrinos is about 600 GeV.    In addition to the total neutrino interaction rates that, for each flavor, are dominated by DIS, we also provide in Table I the expected number of events for several exclusive scattering channels. These include both CC quasi-elastic and NC elastic scatterings (denoted in the table by CCQE and NCEL, respectively), as well as the relevant resonant pion production channels (CCRES and NCRES). We estimate them by convoluting the above neutrino fluxes with the cross sections simulated with GENIE [63,64]. As can be seen, in total approximately 3000 CCQE and CCRES and 1000 NCEL and NCRES events are expected in FLArE-10 during the entire HL-LHC era, and the scattering rate is about 30% larger for FASERν2, and a factor of 7 − 8 larger for FLArE-100. These events are mainly due to interactions of the muon neutrinos, while electron neutrinos are responsible for about 10% of the event rates, and tau neutrinos give subdominant contributions.
As discussed in Ref. [51], the neutrino fluxes predicted by different commonly-used event generators are somewhat different, indicating a flux uncertainty of about a factor of 2. This situation will improve in the coming years, given dedicated theoretical efforts to reduce these uncertainties; see, e.g., Ref. [65]. In addition, measurements of the energy spectra of CC neutrino interactions at FASERν and SND@LHC during LHC Run 3 and later in the FPF neutrino detectors will provide direct measurements of the neutrino fluxes. In the following, we, therefore, assume that the neutrino flux uncertainties are dominated by statistical uncertainties.

C. Signal Modeling
Given our chosen benchmark scenario with m A = 3 m χ , the DM particles originate from the decays of on-shell dark photons produced at the ATLAS IP. We simulate the flux of DM particles through the far-forward detectors with the geometries and locations given in Eqs. (5), (6), and (7), normalizing the number of events to the total integrated luminosity of L = 3 ab −1 anticipated for the HL-LHC era. The dark photons produced in rare π 0 and η decays are obtained by employing the CRMC simulation package [62] and the dedicated EPOS-LHC Monte Carlo tool [66]. In addition, we include dark photon production by dark bremsstrahlung, using the Fermi-Weizsacker-Williams approximation, following the discussion in Refs. [8,36,67].
A rich variety of DM-nuclei scattering processes can be studied with the far-forward detectors. To organize the discussion, in the following, we will divide them into distinct categories in a way similar to neutrino interactions; see Ref. [68] for a review. We first study the case of elastic DM-nucleon scattering, which leads to events with single proton charged tracks in the detector. Next, we consider the exclusive inelastic processes of resonant pion production produced through DM-nucleon interactions. Finally, we consider DM-induced DIS, which is most relevant at highmomentum transfer.

A. Signal
Here we consider elastic DM-nucleon scattering and the associated signature of a single proton track in the detector with no additional visible charged tracks emerging from the interaction vertex. As mentioned above, we will also assume that there is no through-going muon in the detector that could be associated with the DM-induced event. When presenting the results, we will further require that the proton momentum, p p , be above a minimum value defined by the energy threshold of the detector (see Sec. III) and below a maximum value that we chose to maximize the signal to background ratio, S/ √ B. The single proton signature is most directly associated with elastic scatterings of DM off protons, χp → χp. The relevant differential cross section is [31,69] dσ(χp → χp) where E χ is the incoming DM energy, Q 2 = 2m p (E p − m p ) is the squared four-momentum transfer with E p the outgoing proton energy and m p the proton mass, and with τ = Q 2 /(4m 2 p ). The proton form factors can be expressed as where µ p = 2.793, and G D ( As advertised in Sec. II, the scattering cross sections for Majorana fermion and complex scalar DM have the same high-energy limit. This is evident upon inspection of Eqs. (8) and (9), which reveals that the first term proportional to A(Q 2 ) in Eq. (8) is negligible compared to the second term for large E χ . The projected exclusion bounds presented below are therefore valid for both the Majorana fermion and complex scalar DM scenarios. We also note that the integrated cross section becomes independent of the DM energy at large E χ .
Additional signal events could arise from elastic DM scatterings off neutrons, χn → χn, in which the outgoing neutron re-scatters before leaving the nucleus and produces a final-state proton. The relevant cross section for this process can be obtained from Eqs. (8) and (9) by replacing the proton mass and form factors with the quantities appropriate for neutrons [31]. However, because the dark photon mediator couples to electric charge, its coupling to neutrons vanishes in the limit of zero momentum transfer. Therefore, for the models considered here, elastic DM-neutron scattering is considerably suppressed relative to elastic DM-proton scattering. Similarly, inelastic DM scattering followed by the absorption of all charged tracks and neutral pions inside the nucleus, besides a single outgoing proton, contributes subdominantly to the total DM signal event rate. We have verified this using GENIE, under the assumption that the impact of nuclear final-state interactions (FSI) on such particles in DM-induced events can be well approximated by their impact on neutrino events with the same momentum transfer to the nucleus.
In addition to the outgoing proton's energy, its direction can also be observed. Angular cuts were found in Ref. [23] to be useful in separating DM-electron scattering from neutrino-electron scattering, but they are less useful here. In DM-electron scattering, the additional discriminating power was related to the mass hierarchy between the target electron and the incoming DM particles, m e m χ . For the DM-nuclear scattering considered here, however, m χ m p in the parameter space of interest, and so the DM particles behave similarly to essentially massless neutrinos. In the following, we will therefore focus only on the energy cut.
Elastic scatterings χp → χp generally lead to low visible energy depositions due to the strong form factor suppression for large momentum transfers, Q 2 1 GeV 2 . As a result, we will typically set the maximum outgoing proton momentum, p max p , to values below 1 GeV. The DM detection prospects for this signature improve with softer lower limits on the outgoing proton momentum. For this reason, FLArE can be more sensitive than FASERν2 if the FLArE proton kinetic energy threshold, E k,p , can be lowered to 20 MeV, as discussed in Sec. III A. Below, we present in detail the estimated sensitivity reach and background estimates for both types of detectors.
TABLE II. Neutrino-induced background and DM signal events for the single proton signature for several choices of selection cuts on the outgoing proton momentum p p . We assume 14 TeV pp collisions with integrated luminosity 3 ab −1 . The cuts on the minimum proton momentum are dictated by the assumed experimental thresholds, as discussed in Sec. III A. The maximum proton momentum is set to 1 GeV for FASERν2. For FLArE-10 and FLArE-100, we also consider an additional case with p p < 500 MeV. The DM signal corresponds to the benchmark scenario with parameters (m χ , ε) = (100 MeV, 6 × 10 −4 ), m χ = m A /3, and α D = 0.5, and takes into account the efficiency factors (see text).

B. Neutrino-Induced Backgrounds
The dominant neutrino-induced backgrounds to DM-nucleon elastic scattering come from neutral current elastic scatterings (NCEL) of all three neutrino flavors that produce the outgoing proton in the final state, νp → νp. Additional background events can be induced by deep inelastic neutrino scatterings (NCDIS) and resonant pion production processes (NCRES), in which, occasionally, most of the outgoing particles are absorbed in the nucleus due to FSI. We assume below that outgoing electrons and muons can be sufficiently discriminated from protons so that CC neutrino interactions can be neglected in the background discussion.
In Table II, we present the total background event rates obtained with GENIE for FASERν2, FLArE-10, and FLArE-100. In the case of liquid argon detectors, we impose a selection cut on the minimum proton kinetic energy of either E k,p > 20 or 50 MeV, corresponding to the assumed proton detection thresholds discussed in Sec. III A. The latter condition corresponds to a minimum proton momentum of p p 300 MeV, which we also require in the analysis for the emulsion detector. We also cut on the maximum proton momentum, p p < p max p = 1 GeV, and for the more optimistic proton threshold in FLArE, E k,p > 20 MeV, we additionally study a more aggressive upper momentum cut, p max p = 500 MeV. Finally, in each case, we veto on events containing any additional charged tracks or neutral pions emerging from the nucleus, besides the single proton, that have energies above their corresponding detection thresholds; see Sec. III. As can be seen, in the HL-LHC era, the expected number of background events can be roughly 100 events for FLArE-10 and 1000 events for FLArE-100.
The number of background events in FASERν2 is between those in the two liquid argon detectors. The surprisingly large number of expected background events in FASERν2 when compared with FLArE-10, which has a similar mass, is mainly driven by the additional impact of neutrino-induced NCRES events that mimic the single proton signal. The outgoing pions produced in these events often have energies corresponding to the mass difference between the dominant ∆ resonance and the proton, E π ∼ m ∆ -m p ∼ 300 MeV. As a result, such events often lead to pions below the detectability threshold, while the outgoing proton can remain visible. This effect is significantly more pronounced in FASERν2 than in FLArE. A detailed treatment of this background will also depend on the position of the interaction in the tungsten layer, which we leave for future studies with more detailed detector simulations.
For completeness, we also present in Table II the number of DM signal events obtained for a benchmark scenario with m χ = m A /3 = 100 MeV, ε = 6 × 10 −4 (y = 2.2 × 10 −9 ), and α D = 0.5 for three sets of cuts and different detectors. Both in this table and in the subsequent analysis, the number of DM signal events has been additionally rescaled by a finite signal detection efficiency. This is due to the impact of FSI on the outgoing proton that can affect the DM-induced event reconstruction in the detector. We have estimated this efficiency as a function of the momentum of the final-state proton produced in the initial interaction inside the nucleus by studying elastic scatterings of neutrinos with GENIE. The value of the signal efficiency factor that we use in our analysis typically varies between 50% and 70%, and it depends on the energy of the outgoing proton and the analysis type. As can be seen, for FLArE-10 and FLArE-100 with the lower limit E k,p > 20 MeV, the DM signal can yield a 30% to 40% excess over the neutrino background. In contrast, for FASERν2, even though the DM scattering rate is somewhat larger than in FLArE-10, the prospects for probing DM are limited by larger backgrounds.
In the left panel of Fig. 2, we show the signal-to-background ratio S/B as a function of p max p for the FLArE-10 detector. We present results for the above-mentioned benchmark scenario and also for one with (m χ , ε) = (264 MeV, 10 −3 ) (y = 6.2 × 10 −9 ). As evident from Fig. 2, the DM search favors lower values of p max p . This is expected for DM scatterings mediated by the dark photon A , which is much lighter than the Z boson mediating neutrino NC scatterings. For a similar discussion for FLArE and DM-electron scattering, see Ref. [23]. As is apparent from Eq. (8), the lower the A mass, the lower the typical momentum exchange in the χp → χp reaction, which also leads to a lower characteristic momentum of the outgoing proton. In particular, for m A 100 MeV, it would become necessary to require p max p 300 MeV or even lower to obtain S/B ∼ 1. This, however, goes beyond the FLArE and FASERν2 capabilities that we assume in our study. On the other hand, the DM scattering rate can become much higher for increasing mediator mass, in which case a larger momentum exchange is allowed. This can be seen for the case of m A = 3m χ = 792 MeV also shown in the plot. The surprisingly large values of S/B obtained for this benchmark scenario are related to the efficient A production in the proton bremsstrahlung process for m A close to the ρ and ω resonances.
Last but not least, we note that, if systematic uncertainties are negligible relative to statistical uncertainties, the significance of the signal is more closely characterized by S/ √ B than S/B. As p max p increases, the background rate increases, but this increase is milder for √ B than for B, and the dependence on the maximum momentum cut is milder for the ratio S/ √ B than for S/B. For this reason, the projected exclusion bounds shown below are roughly independent of the precise value of p max p .

C. Sensitivity Reach
In the right panel of Fig. 2, we present the expected projected 90% CL exclusion bounds for the three detectors under study. We see that, with just the elastic scattering signature, FLArE-10 will probe most of the thermal relic target for the complex scalar DM model with m χ 100 MeV. For the Majorana fermion DM case, FLArE-10 will only probe the small part of the thermal target region where DM production is enhanced by ω and ρ resonances in the dark photon bremsstrahlung process. The detection prospects could be further improved in the larger FLArE-100 experiment. The expected exclusion bounds for FASERν2 are similar to FLArE-10. We reiterate, however, that, as noted in Sec. III A, this assumes that muon-induced backgrounds can be eliminated for FASERν2 E k,p > 20MeV, p p < 0.5GeV E k,p > 20MeV, p p < 1.0GeV E k,p > 50MeV, p p < 1.0GeV E k,p > 20MeV, p p < 0.5GeV  Table II, and we assume the detectability threshold of E k,p > 20 MeV for the proton kinetic energy. Right: The projected 90% CL exclusion bounds for the elastic scattering signature for FASERν2 with 300 MeV p p 1 GeV (green), FLArE-10 (red), and FLArE-100 (blue) with the proton energy and momentum cuts indicated. Current bounds exclude the gray-shaded region; see Sec. VII for details. The thermal relic targets for the Majorana fermion DM and complex scalar DM models are also shown.

FASERν2.
We also show the impact of different cuts on the proton kinetic energy, E k,p > 50 MeV, and maximum outgoing proton momentum, p max p < 1 GeV. We see that the reach is better in the lowmass region for the lower proton kinetic energy threshold. However, the improved reach mainly corresponds to a region in the parameter space that is already excluded by existing bounds. On the other hand, the expected bound at higher masses is only slightly sensitive to changes of our lower kinetic energy and upper momentum cuts. As a result, the presented sensitivity reach for m χ 100 MeV only mildly depends on the final FLArE capabilities in the considered range of E k,p and p p . When we present combined results for different types of searches in Sec. VII, we will therefore just present results for the cuts E k,p > 20 MeV, p p < 0.5 GeV .

A. Signal
The next signal of interest is χ1π 0 events, in which a single neutral pion is produced through DM-nucleus scattering with no other mesons or charged leptons emerging from the vertex. Such events are produced by DM-induced resonant pion production, χN → χN π 0 , which we model using the BdNMC DM simulation tool [36]. BdNMC accounts for incoherent pion production via excitation of the ∆ resonance, which is expected to be the leading contributor to this process. In addition, χ1π 0 events can also result from DM elastic scatterings off protons followed by FSI. We include this effect in our analysis, although it only mildly affects our final results. When treating the elastic scattering contribution, we assume that the impact of FSI can be modeled using neutrino  [36] that takes into account the dominant pion production via production of the ∆ resonance. We also show the relevant results for neutrino-induced backgrounds from NCRES and NCEL events (brown histogram). This was obtained using the far-forward LHC neutrino energy spectrum and full GENIE [63,64] simulations with further resonances and final-state interactions of hadrons taken into account. Right: The colorful solid lines correspond to the projected 90% CL exclusion bounds in the DMnuclei scattering χ1π 0 signature for FASERν2 (green), FLArE-10 (red), and FLArE-100 (blue). Current bounds and thermal relic targets are as in Fig. 2. interactions, as was discussed in Sec. IV.
In our analysis, we do not differentiate events based on the number of final-state nucleons, including protons, that emerge from the nucleus. This is to mitigate the strong dependence of the number of expected signal events on the assumed FSI model. This inclusive approach is consistent with similar analyses performed by the K2K [70], MicroBooNE [71] and MiniBooNE [72] Collaborations.
The neutral pion in the final state will immediately decay into two photons with momenta typically above the visibility threshold of 30 MeV characteristic for the liquid-argon detectors. In contrast, for FASERν2, the reach will partially be limited by the requirement that photons have an energy of at least 300 MeV to be visible. As discussed above, in the resonant pion production events, we typically observe E π ∼ 300 MeV from the ∆ resonance, which would only be moderately altered by the presence of heavier resonances and FSI. We illustrate this in the left panel of Fig. 3, in which we show the resonant event distribution as a function of the energy of the final-state neutral pion E π 0 for two benchmark DM models with m χ = m A /3 = 10 and 100 MeV, and for neutrinoinduced NCRES background events. The plot has been obtained for the liquid argon detector. As can be seen, in the case of neutrinos, in which the aforementioned effects going beyond the simple ∆ resonance and parton level interactions are taken into account, the resulting distribution is more smeared than for DM. Notably, in both cases, the photons produced in π 0 decays will typically be above 30 MeV.
The characteristic energy of the pions produced through resonant scatterings translates into a relatively weak dependence of the sensitivity reach on the upper energy threshold, which is similar to the elastic DM-nucleon scattering search discussed in Sec. IV. As a result, we will employ a χ1π 0 ν-induced backgrounds Detector E π 0 < 170 MeV 300 MeV 1 GeV 2 GeV  FASERν2  --150  170  FLArE-10  9  90  220  230  FLArE-100  70  740  1750 1850   TABLE III. Neutrino-induced background events in the search for χ1π 0 -type events (see the text for details) as a function of the maximum threshold for the outgoing pion energy. The minimum threshold energy for the outgoing photon is set to 300 MeV and 30 MeV for the emulsion and liquid argon detectors, respectively.
single cut on the maximum pion energy given by E π 0 < 1 GeV. Increasing this limit has a minimal impact on the number of DM-induced resonant pion production events, while it could adversely affect the sensitivity by increasing the number of neutrino-induced backgrounds from DIS events.
Similar to the discussion in Sec. IV, here we also do not discuss the possible impact of the angular cuts on the derived exclusion bounds. We note, however, that the pion angular distribution, as well as the invariant mass reconstruction of the photon pair, could play an important role in further distinguishing such events from neutrino-induced backgrounds producing single electrons in the final state due to the scatterings off electrons or nuclei; see Ref. [38] for a similar discussion for MiniBooNE. Below, for simplicity, we assume that such backgrounds can be rejected in the analysis.

B. Neutrino-Induced Backgrounds
The dominant neutrino-induced backgrounds for the χ1π 0 events are due to NCRES scatterings. We also study subdominant contributions associated with the coherent pion production processes (COHERENT), in which the neutrino scatters off the entire nucleus, and the elastic scatterings NCEL followed by the FSI that generate the outgoing neutral pion. We model all these backgrounds using GENIE. We provide the total expected number of background events for the three detectors in Table III for four choices of the π 0 upper energy threshold: E π 0 < 170 MeV, 300 MeV, 1 GeV, and 2 GeV. As can be seen, increasing the energy threshold above 1 GeV has a very mild impact on the number of background events. We require that the events do not contain any charged pions or other mesons above the visibility thresholds discussed in Sec. III.
Focusing now on the pion energy cut of E π 0 1 GeV, we see that we expect roughly 200 background events in both FASERν2 and FLArE-10, and roughly 2000 such events in FLArE-100. Interestingly, the number of background events is now smaller in FASERν2 than for FLArE-10. This is the opposite effect to the one discussed in Sec. IV, in which increasing the lower energy threshold resulted in a larger number of NCRES events mimicking NCEL ones in the detector. For this reason, we now observe a relatively lower number of NCRES events that will be reconstructed in the emulsion detector as χ1π 0 -like events. As far as liquid argon detectors are concerned, the number of background events in this search is larger, although of a similar order, than for the previously discussed search based on elastic scattering events.

C. Sensitivity Reach
In the right panel of Fig. 3, we present the expected projected 90% CL exclusion bounds based on the χ1π 0 search. As can be seen, the expected bounds are weaker than the ones based on DM elastic scattering shown in Fig. 2. This is primarily due to the smaller scattering cross section. Once we limit the DM signal rate to only NC (A exchange) scatterings off protons with single π 0 production and no charged pions in the final state, the relevant cross section is suppressed relative to elastic scattering by more than an order of magnitude for small mediator masses, m A 100 MeV [36]. The suppression factor becomes smaller, of order a factor of a few, for heavier dark photons. The signal rate is further suppressed by signal efficiencies resulting from FSI and event reconstruction. We estimate them to be of the order of 25% for FLArE and between 10% and 15% for FASERν2. In the latter case, this efficiency also takes into account the aforementioned energy cut of E γ 300 MeV, which is larger than in LArTPC detectors. In the end, we find that the resonant pion signature is less promising than both the electron and single proton signatures.

A. Signal
The last signature that we consider is DM-nuclear scattering at high momentum transfer. Because light DM will be produced with TeV-scale energies in the direction of the FPF, the maximum accessible momentum transfer in nuclear scattering is tens of GeV. Above the QCD scale, deep inelastic scattering leads to a relatively high-energy nuclear recoil, which can subsequently produce multiple charged tracks. In this regime, a partonic treatment is appropriate, and the outgoing hadrons are easily above detector thresholds.
We consider the DIS process χN → χX in the models of Sec. II. The double differential cross section is given by where Q 2 = 2m p E χ xy, x is the parton momentum fraction, y = 1 − E χ /E χ is the fraction of the incoming DM energy transferred to the nucleon in the lab frame, f q is the quark parton distribution function, Q q is the quark electric charge, and As the scattering takes place through a light mediator, it is not surprising that low momentum transfer is favored regardless of the χ spin. Furthermore, the functions B(y) for Majorana fermion and complex scalar DM in Eq. (12) are identical up to O(y 2 ). Because the cross section is dominated by the small y region, then, the DIS signal strength is approximately the same for these two models. This motivates the choice previously mentioned in Sec. II to only show results for the Majorana fermion DM scenario.
To estimate the scattering signal, we convolute these cross sections with the nCTEQ15 parton distribution functions [73] for tungsten and argon nuclei, imposing a minimum cut of Q 2 > 1 GeV 2 . When the parton hadronizes, of course, multiple charged tracks and photons, which yield electromagnetic showers, are produced. We do not simulate this hadronization nor the reconstruction of the hadronic energy and transverse momentum from these objects, though other works have demonstrated the use of track-level information to search for similar signals [74,75]. Instead, we simply take the outgoing parton energy and transverse momentum as proxies for the energy and transverse momentum of the recoiling hadronic system, We expect both E had and p T,had to grow with increasing Q 2 . In principle, there are more detailed kinematic variables involving the visible tracks from the scattered nucleon that could be accessed by doing a full simulation. However, since the hadronic part of each event depends only on the outgoing parton momentum and hadron interaction modeling, we do not anticipate that further kinematic considerations would provide significant additional discriminating power between the signal and neutrino background. The left panel of Fig. 4 shows the two-dimensional distribution of the quantities in Eq. (13) for DIS in one of our benchmark DM scenarios at FLArE-10. The distribution is qualitatively similar at FASERν2. The signal events are clustered at lower energies and transverse momenta than the background, consistent with the preference for low momentum transfer in light DM scattering. Despite the preference for low momentum transfer, there is still a significant number of events with energetic nuclear recoils. We see the most events at E had of several GeV and low p T,had , and expect that such events would have multiple tracks emerging from a vertex with no incoming track. A more detailed study of the detection efficiency, including the effects of hadronization and FSI, would be interesting. For instance, the efficiency would depend on the number of tracks and hence the hadron multiplicity, which tends to grow with the center-of-mass energy W of the recoiling hadronic system. W is related to the momentum transfer and partonic momentum fraction through W 2 = m 2 p + Q 2 (1 − x)/x. The EMC experiment measured the charged hadron multiplicity in muon DIS, finding that several charged tracks were typical for W > 4 GeV [76]. We have checked that a cut of W > 2 GeV, which would avoid the resonant scattering region with fewer tracks, does not change our results significantly. In addition, as our signal is clustered at values of p T,had /E had corresponding to angles of several degrees, it would be useful to examine technologies for measuring multiple hadronic tracks in the forward direction in liquid argon for FLArE. While there can be difficulties measuring such tracks using wire planes if the planes are oriented parallel to the track direction, the patterns of charge deposition can be used to obtain three-dimensional information [77], as has been demonstrated by MicroBooNE for neutrino event identification [78] and cosmic ray rejection [79].

B. Neutrino-Induced Backgrounds
The main background to DM DIS is neutrino scattering. NC neutrino scattering would produce a nuclear recoil with significant energy carried away by the outgoing neutrino, just as in our signal. CC neutrino scattering, by contrast, would result in a high-energy outgoing lepton. We assume that the detector would have sufficient efficiency that the neutrino CC backgrounds could be rendered very small.
There are also backgrounds from muon interactions, which can be eliminated by requiring that there is no charged track leading into the vertex [13]. Muon interactions can also produce neutral hadrons, which travel for distances on the order of 10 cm before scattering. These neutral hadron events can mimic the signal. Although neutral hadron backgrounds are problematic in a pure emulsion detector [13,75], as mentioned in Sec. III A, we assume that an active muon veto will remove these events at FASERν2 or FLArE [23]. By using timing to remove a small area around each muon interaction, we expect that neutral hadron scattering could be reduced to negligible levels without significant impact on the signal.
The differential NC neutrino scattering cross section at high energy is [80] dσ(νN → νX) dx dy in terms of the partonic momentum fraction x and the fractional neutrino energy loss y = 1 − E ν /E ν = E had /E ν . The momentum transfer is Q 2 = 2m p E ν xy. Here, g q,L , g q,R = T 3 − Q sin 2 θ W are the NC couplings of the quarks. For anti-neutrinos, the cross section is As the momentum transfer Q 2 is generally small compared to m 2 Z , the neutrino scattering cross sections are proportional to the CM energy or, equivalently, the energy of the incoming neutrino.
The typical Q 2 is perhaps the most striking difference between light DM DIS and neutrino NC scattering. In principle, the momentum transfer 2E χ m N in DM scattering can be as high as tens of GeV. However, for scattering through a light mediator, smaller momentum transfers are typically preferred, as the scattering cross section goes as 1/Q 4 in the limit of vanishing mediator mass. On the other hand, neutrino scattering proceeds through the Z, which is heavy compared to the typical momentum transfer. Consequently, the neutrino NC scattering cross section grows linearly with the partonic CM energy √ŝ . We proceed to investigate the kinematics further to discriminate between signal and background, showing the hadronic energy and transverse momentum for the neutrino background in the right panel of Fig. 4 and p T,had : Strong cuts: 1 GeV < E had < 15 GeV , 1 GeV < p T,had < 1.5 GeV Loose cuts: 1 GeV < E had < 30 GeV , 1 GeV < p T,had < 2.0 GeV .
The effects of these cuts on the background and signal are shown in Table IV. We see that the background can be reduced by over an order of magnitude while keeping 1/4 to 1/2 of the DM DIS signal.

C. Sensitivity Reach
Having examined the kinematics of the signal and background events, we present the expected projected 90% CL exclusion bounds for DM DIS searches at FASERν2 and FLArE in Fig. 5. Considering statistical uncertainties only, the former set of cuts in Eq. (16) yields the strongest projected exclusions. The figure shows the reach of the different detectors, as well as the effects of the hadronic energy and transverse momentum cuts in the case of FLArE-10. In contrast to the lower energy signatures in Secs. IV and V, the typical deposited energy is well above the thresholds for both emulsion and liquid argon detectors. The relative performances of FASERν2 and FLArE thus depend mostly on the detector mass and geometry, as well as on their background rejection and event identification capabilities. Here, we focus on the former, while assuming 100% signal detection efficiency for both types of experiments. Of the two 10-tonne detectors, the more compact FASERν2 provides better sensitivity to light DM scattering because it has more mass at large rapidity where the DM flux is higher. In addition, the numbers of events for FLArE-100 in Table IV do not scale fully with the detector mass, when compared to its 10-tonne analog. Similar effects were observed for DM-electron scattering [23].
As discussed above, the DIS limits are very similar for the Majorana fermion and complex scalar DM models, and we have used the former to draw the projected exclusion lines. To guide the interpretation of the limits, we also show the thermal relic targets in each of these scenarios, assuming standard thermal cosmology. We see that DIS searches can probe dark photon scenarios yielding the correct thermal relic density for DM masses above approximately 200 MeV. The expected sensitivity reach can then also partially cover the resonance region, in which the intermediate dark gauge boson in DM annihilations mixes with the SM vector mesons ρ and ω, i.e., 2m χ ≈ m ρ,ω , especially for complex scalar DM. By contrast, the reach of DM DIS is relatively limited at low masses. This is because the growth of the DIS cross section at small mediator masses is limited by our minimum momentum transfer cut of 1 GeV. Nevertheless, DM DIS searches at FPF detectors offer the potential to probe dark photon scenarios that are viable from the standpoint of thermal cosmology and that are otherwise unconstrained. Finally, we note from Table IV that with the full HL-LHC dataset, there will be thousands of background events even with kinematic cuts. It will thus be important to reduce uncertainties from systematics such as the neutrino flux and signal/background modeling, which we have not considered here, in a full experimental analysis. We assume that they can be suppressed so that the analysis will be dominated by statistical uncertainties. For instance, as has been suggested previously [23], measuring the neutrino flux at other detectors or in other kinematic regions could help constrain the background normalization. If statistical uncertainties dominate, then since the number of signal events scales with y 2 , the limit on y associated with a fixed significance S/ √ B improves as L −1/4 . The impact of this mild dependence is that new parameter space can be probed with a relatively small amount of data. We will consider the effect of luminosity on the reach more completely in the next section, where we combine the results of this section and the previous two to obtain the overall FPF reach in searches for light DM-nuclear scattering.

VII. COMBINED SENSITIVITY REACH
In this section, we combine all of our previous results on DM-nucleus scattering processes, including elastic scattering (Sec. IV), resonant pion production (Sec. V), and DIS (Sec. VI), as well as the results previously obtained [23] for the DM search based on scatterings off the electrons. These are shown for FLArE-10 in Fig. 6. In general, since the scattering cross sections grow for small mediator mass and we have taken a fixed mass ratio m A /m χ , the limits are strongest at low m χ . The flattened sensitivity at the left of the plot arises from the minimum momentum transfer for each process considered. For elastic scattering and resonant production, these come from experimental considerations on the visibility of the outgoing proton or pion. We see that the low thresholds of liquid argon detectors allow for the ability to probe new parameter space at m χ 200 MeV. For DIS, the Q 2 cutoff to ensure the validity of our partonic treatment limits the sensitivity at small masses, but the inherently harder nature of DIS can lead to stronger bounds at higher DM mass. Figure 6 also shows that elastic scattering and DIS are the most sensitive nuclear scattering probes at low and high masses, respectively. Resonant pion production is never the strongest signature in this model. The sensitivity reach from DM-electron scattering, derived previously in Ref. [23], is also shown, and can be seen to be competitive with the best nuclear signatures at moderate and high masses, and even stronger at low masses.
In Fig. 6, we also show the thermal relic targets for Majorana fermion and complex scalar DM, as well as current and projected results from other experiments. Existing bounds from null results are shown as the gray shaded region. These include results from BaBar [82], MiniBooNE [38], and NA64 [83], as well as recasts of searches at BEBC [84], CHARM-II [85], E137 [86,87], LSND [29], and NOνA [88], as derived by the authors of Refs. [81,89]. Projected sensitivities of future experiments are shown in the dashed and dotted colored contours. We also note that future short baseline neutrino experiments such as ICARUS could also be sensitive to DM scattering [81]. The brown contours are projected sensitivities from searches for DM that is produced at a collider or beam dump and then subsequently scatters in a downstream detector, a signature similar to what we have considered in this work. These include BDX [90], SND@LHC [53], and SND@SHiP [91]. The red contours are projected sensitivities of future missing momentum-type searches, including NA64 [92], LDMX [47,93], and Belle-II [94]. Last, the green contour shows the projected sensitivity of SuperCDMS to the Majorana fermion DM model [42,47,93]. The region probed by SuperCDMS is at higher masses than are probed by FLArE-10. For the complex scalar DM model, direct detection limits can be more constraining, but they can also be evaded by the introduction of a small mass splitting between the DM states so that the scattering is inelastic. Figure 7 then shows the best limits from each of the detectors in Sec. III. As for FLArE-10 in Fig. 6, the best limits arise from electron scattering and nuclear DIS in the low and high mass ranges, respectively. Both FASERν2 and FLArE-10 will probe the relic target for the complex scalar DM model for DM masses between several MeV and a few hundred MeV. FLArE-100 could provide a similar reach for the Majorana fermion DM model. Altogether, the detectors we have studied are able to probe a wide swath of the cosmologically-favored parameter space for both the Majorana fermion and complex scalar DM models.
Finally, to estimate the time scales on which forward LHC searches could start to achieve new sensitivity to light DM, we show the 90% projected exclusion bounds at FLArE-10 for a selection of integrated luminosities in Fig. 8. Again, the best limits from all processes (elastic proton scattering, resonant pion production, DIS, and electron scattering) have been used. With around 30 fb −1 of data, these searches can begin to test thermal DM scenarios that are thus far unconstrained. In addition, the 5σ discovery reach as a function of m χ is a factor of approximately 1.6 in y above the projected exclusion bounds. As a result, DM can be discovered at the 5σ level with 3000 fb −1 for DM masses of 3 -10 MeV and 50 -300 MeV.

VIII. CONCLUSIONS
The search for terrestrial DM production is a major component of the physics programs of particle accelerator and collider facilities. This avenue is especially useful in the case of sub-GeV DM, where traditional direct detection experiments lose sensitivity. Such light DM at the LHC would be dominantly produced at high rapidities beyond the reach of the general-purpose detectors, motivating dedicated experiments in the far-forward direction. In this work, we have studied potential DM scattering signals in forward detectors at the HL-LHC, as would be possible at the FPF. Our focus has been on interactions between DM and nuclei, complementing previous work on DM-electron scattering.
We have considered detectors based on both emulsion and liquid argon technology. With the use of timing information, it would be possible to reject muon-induced backgrounds, including those from neutral hadron interactions. We thus expect that the dominant backgrounds will be from neutrino scattering. Indeed, the scattering processes that we have considered are analogous to SM processes with neutrinos: elastic scattering, resonant pion production, and DIS. For each of these processes, we have estimated the DM signal and neutrino background, investigating the differences due to kinematics and incorporating the effects of nuclear FSI as appropriate. We find that for DM scattering through a light mediator, it is possible to mitigate neutrino backgrounds with kinematic cuts favoring events with low momentum transfer. This strategy is effective because the heavier weak gauge bosons cause neutrino backgrounds to prefer high Q 2 scattering. Similar considerations apply to other signatures, and it would be interesting to study whether additional sensitivity could be obtained with other processes. These include coherent scattering, coherent pion production, and multiple meson production.
Looking at benchmark models with light DM interacting through the minimal dark photon portal, we find new sensitivity in the MeV to GeV mass range. With either complex scalar or Majorana fermion DM, the searches here would test regions of parameter space in which the observed relic density is obtained through thermal freeze-out. As the characteristic energies of the processes that we have studied are different, they have complementary sensitivities. When these searches are combined with those for DM-electron scattering, FASERν2 and FLArE-10 could cover the relic target for complex scalar DM for DM masses between several MeV and several hundred MeV. FLArE-100 would provide sensitivity to the relic target in a similar mass range for Majorana DM, which is not probed by current experiments. All of these experiments cover much of the parameter space in which the thermal relic density does not overclose the Universe, and they have the potential to provide direct evidence for DM interactions, in contrast to missing momentum-based searches at accelerator and beam dump facilities. Notably, currently unconstrained regions of parameter space can start to be probed with even the first O(30) fb −1 of integrated luminosity at the HL-LHC.
The forward region of the LHC offers exciting possibilities to study physics within and beyond the Standard Model. The FPF would extend the reach of the LHC, providing qualitatively new discovery potential in well-motivated theories of light dark sectors. In addition to electron scattering, a suite of nuclear scattering searches at the FPF detectors can be performed to improve our understanding of the nature of DM. In searching for DM and beyond, further exploration of the unique environment provided by the far-forward region at the LHC is warranted to fully leverage collider probes of new physics.