Hadrophilic Dark Sectors at the Forward Physics Facility

Models with light dark sector and dark matter particles motivate qualitatively new collider searches. Here we carry out a comprehensive study of hadrophilic models with U(1)$_B$ and U(1)$_{B-3L_{\tau}}$ gauge bosons coupled to light dark matter. The new mediator particles in these models couple to quarks, but have suppressed couplings to leptons, providing a useful foil to the well-studied dark photon models. We consider current bounds from accelerator and collider searches, rare anomaly-induced decays, neutrino non-standard interactions, and dark matter direct detection. Despite the many existing constraints, these models predict a range of new signatures that can be seen in current and near future experiments, including dark gauge boson decays to the hadronic final states $\pi^+ \pi^- \pi^0$, $\pi^0 \gamma$, $K^+ K^-$, and $K_S K_L$ in FASER at LHC Run 3, enhancements of $\nu_{\tau}$ scattering rates in far-forward neutrino detectors, and thermal dark matter scattering in FLArE in the HL-LHC era. These models therefore motivate an array of different experiments in the far-forward region at the LHC, as could be accommodated in the proposed Forward Physics Facility.


I. INTRODUCTION
Searches for new particles and dark matter (DM) are primary physics drivers at the Large Hadron Collider (LHC). Traditional searches for the classic missing p T signature at the LHC main detectors have sensitively searched for particles with weak-scale masses and O(1) couplings to Standard Model (SM) particles, but are less effective for light and weakly coupled new particles, including long-lived particles (LLPs) and DM. Recently it has been appreciated that new experiments in the farforward region at the LHC can provide a powerful probe of new light particles. These experiments exploit the large forward flux of pions and other SM particles, which, if they decay to new light particles, can create a large forward flux of LLPs and DM. Light new physics species can also be produced in the far-forward region of the LHC in other types of interactions, including proton-proton bremsstrahlung and the Drell-Yan process. The recent detection of TeV neutrino candidates in the forward re- * batell@pitt.edu, jlf@uci.edu, mfieg@uci.edu, aismail3@okstate.edu, felix.kling@desy.de, rmammen@okstate.edu, strojanowski@camk.edu.pl gion [1] also opens a new window on neutrinos at colliders, which may be used to probe both SM and beyond the SM (BSM) phenomena [2][3][4].
In evaluating any proposal for new physics at the MeV to GeV mass scale, one must carefully consider all of the existing constraints from particle and nuclear experiments carried out over the last 60 years. To do this requires a model framework. The dark photon model has been discussed at length in the literature. It is theoretically attractive and contains within it phenomenologically-viable benchmark scenarios of light thermal DM. Of particular relevance for this study, previous studies in the dark photon framework have established the potential for forward experiments to detect both LLPs [5,6] and light thermal DM [7,8]. At the same time, the experimental signatures of a given dark sector model are, to a large extent, determined by the interactions of the mediator with the SM. To more fully evaluate the physics potential of proposed experiments, then, a variety of phenomenologically distinct mediators must be examined. Since the LHC is a pp collider, it is natural to consider mediators with hadrophilic couplings, i.e., sizable couplings to quarks, but suppressed couplings to leptons. Although such models are challenging to test at electron facilities (e.g., Belle-II [9], NA64 [10], 2 LDMX [11], and SENSEI [12]), one might suspect that they can be sensitively probed at proton facilities, such as the LHC.
In this work we study the prospects for probing two dark sector models with hadrophilic vector boson mediators. The first model is based on a gauged U(1) B baryon number symmetry (see, e.g., Refs. [13][14][15][16][17]). This model is perhaps the first example of a hadrophilic model one might consider, since it has sizable couplings to quarks and (loop-)suppressed couplings to all leptons. The model suffers from gauge anomalies, however, which potentially lead to stringent constraints from rare FCNC and Z boson decays [18,19]. We will evaluate the prospects for discovering new physics in this model, carefully respecting all anomaly constraints, as well as those from other experimental searches. We note that anomaly-free extensions of the SM with a local U(1) B symmetry and DM have been constructed in Refs. [20][21][22][23], which focus on the case of new particle masses above the weak scale.
As a second example we consider a model with a U(1) B−3Lτ vector boson mediator. (In the rest of this paper, we will use the modest abbreviation of U(1) B−3τ for this symmetry.) With the addition of a right-handed neutrino, this symmetry is anomaly free and therefore evades the most stringent rare decay constraints present in the U(1) B model. This model is also hadrophilic, in the sense that couplings to electrons, muons, and their accompanying neutrinos are suppressed. However, the presence of τ and ν τ couplings brings with it both additional constraints from neutrino nonstandard interactions (NSI), and also new opportunities for signals involving the 3rd generation leptons. A goal of this study is to incorporate all these new constraints and see what discovery prospects remain.
We will consider both current and proposed farforward experiments. In the last two years, the magnetic spectrometer and tracking detector FASER [24], and the two emulsion detectors FASERν [25] and SND@LHC [26] have been approved. FASER has been fully constructed, and all three are expected to begin taking data when Run 3 starts in 2022. For the High Luminosity LHC (HL-LHC) era, detectors under consideration include upgrades of these detectors (FASER2, FASERν2, and Advanced SND), as well as the Forward Liquid Argon Experiment (FLArE) [7]. 1 A new facility, the Forward Physics Facility (FPF) [4,29], has been proposed to accommodate these experiments.
Remarkably, we will find that all of these detectors have discovery prospects for the hadrophilic models we consider. The possible signals include DM deep inelastic scattering (DIS) and elastic scattering, enhanced predictions for neutrino neutral current (NC) scattering, an excess of tau neutrinos in the forward region, and the visible decay of the dark mediators into SM final states. Notably, the visible decays include final states, such as π + π − π 0 , π 0 γ, K + K − , and K S K L , that could conceivably appear in FASER at LHC Run 3; such states are inaccessible at FASER in dark photon models. The signals are diverse and require a similarly diverse set of experiments to find them, and when combined, the experiments probe parameter space even beyond the DM thermal targets. These models therefore add to the broad physics portfolio of the FPF, complementing other studies of long-lived particle searches, collider-produced TeVenergy neutrinos, new probes of QCD, and high-energy astroparticle physics [4].
The paper is organized as follows. In Sec. II we introduce the two hadrophilic dark sector models based on the U(1) B and U(1) B−3τ gauge symmetries and discuss the production and decays of the vector boson mediator, the DM thermal relic abundance, and the existing constraints for each model. Next, we present our assumptions regarding the performance of FASER, FASER2, SND@LHC, FASERν2, and FLArE in Sec. III. In Sec. IV we outline our methodology for estimating the sensitivity of these far-forward detectors to the new physics signatures predicted in these hadrophilic models. Our main results are contained in Sec. V, and our conclusions and outlook are presented in Sec. VI.

A. Models
With the motivation outlined in Sec. I, we begin in this section by describing the two representative hadrophilic dark sector models based on the anomalous U(1) B and anomaly-free U(1) B−3τ gauge symmetries. 2 Since the new gauge group is Abelian, the new vector gauge boson generically mixes with the SM photon through a kinetic mixing term F µν V µν , where F µν and V µν are the field strengths of the SM photon and new gauge boson, respectively. In the physical mass basis, the Lagrangian of the vector boson mediator V µ is where m V is the vector boson mass, J µ SM is a current composed of SM fields, and J µ χ is the current for the dark matter particle χ.
The SM current is EM are the baryon number and electromagnetic currents, respectively, ε is the kinetic mixing parameter, and x = 0 (1) for the U(1) B (U(1) B−3τ ) model.
To specify J µ χ , we must choose the DM candidate χ. We will study both complex scalar DM and Majorana fermion DM in this work, with Lagrangians where m χ is the DM mass. The associated currents, J µ χ in Eq. (1), are where Q χ is the charge of the DM under the new gauge symmetry. As we will discuss below, both complex scalar and Majorana fermion DM exhibit velocity-suppressed P -wave annihilation to SM final states, implying that bounds from precision measurements of the cosmic microwave background anisotropies [31,32] are easily satisfied in these models. Furthermore, Majorana DM features momentum-dependent scattering in the nonrelativistic regime, making it challenging to probe with DM direct detection experiments. This is not the case for complex scalar DM, and, as we will see, direct detection experiments place strong constraints on such DM for masses above the GeV scale. However, it is important to note that these constraints can also be evaded in a straightforward way by introducing a small mass splitting, which renders the scattering transition inelastic [33][34][35].
The full parameter space of these models is, then, specified by 5 parameters: To reduce the parameter space, as is commonly done in the literature, we will assume a kinetic mixing parameter of typical one-loop size, This is the parametric size of the kinetic mixing generated by loops of SM particles charged under both electromagnetism and the new gauge symmetry. The kinetic mixing depends, in general, on the details of the UV physics and therefore cannot be determined unambiguously, but we neglect such effects here; see also Ref. [36] for further discussion of this issue. Throughout our study we will also adopt another common convention, so that DM annihilation proceeds to SM particles through a virtual s-channel vector boson mediator. Given the assumptions of Eqs. (6) and (7), the resulting parameter space may be specified by the three parameters We will present our results in the (m V , g V ) plane with various choices for Q χ . Since the new symmetries are Abelian, the charge Q χ may be any real number. When presenting our results below, we will consider two choices for coupling hierarchies. As a first scenario, we will consider DM and SM particles to have comparable interaction strengths with the vector boson mediator, fixing In the B − 3τ model, we have fixed the DM charge to be opposite that of the ν τ , Q χ = −Q τ . As a second, qualitatively distinct, scenario, we consider the case in which the DM coupling to the vector boson mediator has a fixed value, Given that we will be considering vector boson mediators with weak couplings to the SM, that is, values of g V ∼ 10 −8 − 10 −2 , Eq. (10) implies very large DM charges Q χ . This may appear unnatural, but there is nothing wrong in principle, since the expansion parameter α χ remains perturbative. Ideas for achieving such large coupling hierarchies for two U(1) gauge symmetries have been presented in Ref. [37]. Finally, although we do not consider them in this work, viable models of hadrophilic scalar mediators can also be constructed; see, e.g., Refs. [38][39][40]. However, for the incident DM energies in the TeV range relevant for FPF experiments, scalar-mediated DM-nuclear scattering rates are typically suppressed by several orders of magnitude in comparison to vector boson-mediated scattering rates; see also Ref. [37] for a comparison of vector boson-and scalar-mediated DM scattering in the ultra-relativistic regime. For this reason, scalar-mediated DM scattering can be better probed by low-and medium-energy experiments [39]. On the other hand, experiments such as FASER and FASER2 can have powerful sensitivity to visible decays of the long-lived scalar mediator in these models, as has been demonstrated in Ref. [41].

B. Production and Decay of the Vector Boson Mediator
In our simulations, we model the production of light dark vector bosons in the far-forward region of the LHC by employing the FORESEE package [41]. We thereby include dark vector boson production by light meson decays, proton bremsstrahlung, 3 and the Drell-Yan process. We observe that typically the production of dark vector bosons in light meson decays dominates if kinematically allowed. For the dark vector bosons heavier than the η meson, the most important production mode is due to bremsstrahlung, while the Drell-Yan process starts to dominate for m V > 1.5 GeV.
We then consider various decay final states of the dark vector bosons. In particular, the partial decay width for V → χχ * is where κ = 1 and 2 for complex scalar and Majorana DM, respectively. The partial decay width into hadrons and other SM particles is taken from the DarkCast package [45], which used data-driven methods to estimate the hadronic width. An alternative description has also recently been implemented in Herwig 7; see Ref. [46]. In Fig. 1, we present the corresponding decay branching fractions for both of the models assuming the Q χ charge as in Eq. (9). In the case of the U(1) B model, LLP decays into lepton pairs are always subdominant, since they appear only at the loop level through the vector boson mixing with the photon. In contrast, the invisible branching fraction of V → χχ * is close to unity for light vector boson masses up to the ω-resonance region, m V ≈ m ω 782 MeV. This leads to an intense flux of DM particles, which can be detected via DM scatterings, as we will discuss in Secs. IV A and IV B. For heavier dark vector bosons, decays into light hadrons start to play an important role and can lead to additional signatures in the detectors, as we will see in Sec. IV E.
For the B − 3τ model with the dark charge set to Q χ = 3, we obtain BR(V → χχ * ) ∼ (10 − 20)% up to the tau threshold, above which V → τ + τ − decays become kinematically allowed. The remaining decay rate for lighter dark vector bosons is dominantly into tau neutrinos, V → ν τ ν τ . As will be discussed in Sec. IV D, this can contribute to the total ν τ flux measured at the FPF. The decays into hadrons also become important for certain values of m V , especially around the ω-and φ-resonance regions.

C. Thermal Relic Abundance
Thermal targets, that is, the regions of parameter space where DM annihilates in the early Universe through thermal freezeout to the correct relic density, provide an important standard by which to judge the sensitivity of collider searches. These have been determined in the U(1) B model with fixed α χ = 0.5 in Ref. [37]. Here we determine, for the first time, the thermal targets for the U(1) B model with fixed Q χ and for the U(1) B−3τ model described above.
The dark matter annihilation cross section can be written in the standard resonance form, where β χ (s) = (1 − 4m 2 χ /s) 1/2 , s V = 1, s χ = 0, and Γ χχ * (s) and Γ SM (s) are the partial decay widths for V decaying into dark matter and SM particles, respectively, with the replacement m V → √ s. The thermally-averaged cross section is, then, [47] where v is the relative velocity of the annihilating dark matter particles, and K i is the modified Bessel function of order i. To determine the thermal target regions of parameter space, we require which reproduces the observed DM relic abundance for the masses we consider [48]. The thermal targets are presented below in Figs. 2 and 3. Their shapes can be understood as follows. In the U(1) B−3τ models, annihilation to tau neutrinos is allowed throughout the m V range. The thermally-averaged cross section has the parametric dependence and so in the (log m V , log g V ) plane, the thermal targets have slope 1 for the models with fixed α χ shown in Fig. 2, and slope 1/2 for the models with fixed Q χ shown in Fig. 3. The discrepancy between the complex scalar and Majorana fermion cases results from the fact that in the complex scalar case, there are both DM and anti-DM particles, whereas in the Majorana case, DM is its own anti-particle, which impacts the annihilation rate through the parameter κ's appearance in Eqs. (12) and (13). The "heavy hadrons" contour includes charm and bottom hadrons, and the red contours correspond to all other hadrons. Among them, we explicitly show the dominant branching fractions into π 0 π + π − , π 0 γ, and kaon pairs KK = K + K − + KSKL. Here we assume loop-induced couplings of the bosons to charged leptons of the first two generations of size g = gV (e/4π) 2 . The relevant contour for boson decays into e + e − or µ + µ − , shown in the left panel, has been multiplied by a factor of 1000 for visibility. The DM is taken to be a scalar, with the decay width given in Eq. (11).
For the U(1) B models, the thermal target slopes are similar to those for the U(1) B−3τ models for m V 1 GeV. The required couplings g V are greater because the annihilation to tau neutrinos is absent. As m V drops below 1 GeV, the cross section to hadrons decreases rapidly, and without a large leptonic annihilation channel, the required g V increases rapidly to maintain a fixed σ ann v . This continues until m V drops below m π , at which point all hadronic channels shut off, and only the loop-suppressed annihilation to light leptons is allowed. The curve moves further up for masses m V /3 = m χ < m e where only the high velocity tail of the thermal DM population can annihilate into electrons, which needs to be compensated by a larger coupling. However, even though only a small fraction of DM can annihilate into electrons, this is still more efficient than the annihilation into 3 photons. The latter process, χχ → 3γ [49,50], was found to be negligible for our study.
The resonance structure seen in all cases arises from resonant mixing of the dark gauge boson V with the SM vector mesons ρ, ω, and φ. In the case of DM annihilation, these resonances occur at masses 2m V /3 = 2m χ = m ρ,φ,ω , whereas for V production, these resonances occur at m V = m ρ,φ,ω .

D. Existing Constraints
Light hadrophilic mediators have a rich phenomenology, giving rise to constraints from previous searches, as well as search opportunities at FPF experiments. Below, we summarize the various laboratory experimental constraints on light hadrophilic gauge bosons following the discussion of Ref. [51]. The resulting limits are shown in Figs. 2 and 3 as dark gray shaded regions.
Invisible Mediator Decays: The focus of this study are hadrophilic mediators with a sizable branching fraction into dark matter. This decay leads to missing energy signatures which have been searched for by various experiments. The most sensitive constraints have been obtained by the search for the decay π 0 → γV at NA62 [52] and LESB [53]; the search for the decay π 0 , η, η → γV at Crystal Barrel [54]; the search for the decay K + → π + V at E949 [55] as discussed in Ref. [15,56]; the search for the mixing induced invisible decays of the J/Ψ by BES [57] and the Υ by BaBar [58] as discussed in Ref. [15]; and the monojet search pp → V + jet at CDF [59] as discussed in Ref. [60].
Visible Mediator Decays: If the couplings of the hadrophilic mediator to the SM and dark sector have similar size, decays into visible final states are possible. If the coupling is sufficiently large, the decays of the mediator occur promptly in the detector and can be searched for via a bump hunt. Bounds have been obtained by the search for the decay η → V γ → π 0 γγ at GAMS-2000 [61] and the search for the nonelectromagnetic contribution to the decay Υ(1S) → jj by ARGUS [62], as discussed in Ref. [63]. In addition, there are bounds from searches for displaced decays of LLPs from NuCAL [64].
DM and Neutrino Scattering: The hadrophilic mediator is copiously produced in beam dump experiments. The decay V → χχ * then leads to a dark matter beam. The MiniBooNE collaboration has searched for the scattering of χ in their downstream neutrino detector [65,66]. Recently, even stronger bounds on coherent scatterings of leptophobic DM have been obtained with the Coherent CAPTAIN-Mills (CCM) liquid argon (LAr) detector [67]. Similarly, the decay V → ν τ ν τ leads to an increased tau neutrino flux, which can be 6 constrained using measurements from DONuT [68], as discussed in Ref. [51].
Indirect Probes: A hadrophilic mediator can also be constrained indirectly through its contribution to the low-energy neutron-lead scattering cross section [69], as discussed in Ref. [70]. Additionally, a new gauge boson with couplings to tau leptons can be constrained by the measurement of the Z → τ τ decay width at LEP [71], as discussed in Ref. [72].
In addition, there are other constraints that are somewhat more model dependent. These are the anomaly constraints and the constraints from neutrino NSIs, which are shown as light gray shaded regions in Figs. 2 and 3, and which we now describe: Anomaly Constraints: As mentioned above, the dark vector boson in the U(1) B model couples to a nonconserved SM current. Invisible decays of such a vector boson are then constrained by enhanced bounds from missing energy searches in rare Z decays and flavorchanging meson decays K → πV and B → KV . We implement them following Refs. [18,19], assuming that anomalies associated with the new gauge group are canceled by heavy fermions that do not receive masses from electroweak symmetry breaking. If these anomalies were canceled by fermions with Yukawa couplings to the Higgs, the invisible decay constraints would not apply, but there would be severe LHC constraints on the additional fermions.
Neutrino NSI: For the U(1) B−3τ model, additional constraints arise from studying neutrino oscillations, both in vacuum and in the matter background of the Sun and Earth. These have been precisely measured by a variety of neutrino experiments. A global fit to these neutrino oscillations measurements simultaneously constrains the oscillation parameters and NSI between neutrinos and matter. We present these bounds following Ref. [73]. 4 We note, however, that these constraints are model dependent and could be weakened in the presence of additional new physics.
Direct Detection: Further bounds on hadrophilic DM can arise from direct detection (DD) searches [15]. These, however, depend sensitively on the detailed structure of the DM interaction and do not apply to Majorana DM and to inelastic scalar DM if the mass gap between the dark species is large enough to suppress upscatterings of non-relativistic DM particles. We stress this in the following when presenting the current DD bounds on spin-independent DM-nuclei scattering from the CRESST-III [76], DarkSide-50 [77], and Xenon 1T [78,79] experiments. We show these bounds assuming that Ω χ h 2 0.12 [80] in the entire reach plot and that a non-standard cosmological scenario affects the DM relic density for points in the parameter space away from the thermal target lines.
Cosmology & Astrophysics: Further indirect probes arise from possible contributions of light dark vector bosons to the number of relativistic degrees of freedom in the early Universe, ∆N eff . We present them below following Refs. [81,82]. Additional bounds could arise from an enhanced supernova cooling rate of SN1987A, as discussed, for example, in Refs. [83][84][85][86][87][88][89]. Such constraints typically probe very small couplings outside the regions of interest for this study. In addition, they are also dependent on a number of astrophysical assumptions, which may weaken the constraints or possibly even remove them altogether; see, e.g., Ref. [90]. In the following, we do not show these bounds explicitly in our sensitivity reach plots, as a detailed study for the models considered here is beyond the scope of our analysis.

III. DETECTORS
We perform our analysis for the on-axis far-forward detectors that will operate either during LHC Run 3 or the HL-LHC era. In the latter case, we focus on the proposed FPF, which begins at a distance L = 620 m away from the ATLAS Interaction Point [4]. In particular, we study the expected future sensitivity of the 10-tonne emulsion detector FASERν2, a proposed successor to the FASERν experiment that will take data during LHC Run 3 [2,25], as well as the 10-and 100-tonne fiducial mass liquid-argon time projection chamber (LArTPC) detectors FLArE-10 and FLArE-100 [7]. The relevant detector geometries are where ∆ is the length of the detector, and S T denotes its transverse size.
Both types of detectors have excellent capabilities to reconstruct the low-energy nuclear scattering signals created by both neutrinos and hadrophilic DM, and also to disentangle DM-induced events from the more energetic neutrino scatterings. For these searches, however, it is also important that they be able to reject backgrounds induced by high-energy muons that pass through the facility and interact with the surrounding rock and infrastructure. To veto these muons, it is highly beneficial to collect time information about the events. In the case of FASERν2, this would likely require interleaving the emulsion layers with additional electronic detectors. For FLArE, on the other hand, the required time resolution can be more easily obtained by employing an additional light collection system; see Ref. [7] for further discussion.
Throughout this paper, we use the neutrino fluxes for the FPF as presented in Ref. [8]. These fluxes were obtained using the event generator SIBYLL 2.3d [91][92][93][94] as implemented in CRMC [95] to simulate the primary collision, and the fast neutrino flux simulation presented in Ref. [96] to model the propagation and decay of long-lived SM hadrons in the forward LHC infrastructure.
In addition to the aforementioned scattering detectors, we will also present sensitivities for the LLP signature of the vector boson mediator decaying to visible SM final states. To this end, we will focus on FASER [24,97] and FASER2, cylindrical detectors with length ∆ and radius R, where FASER will take data during LHC Run 3 and will be positioned in the far-forward region at a distance L = 480 m away from the ATLAS IP. For FASER2 we assume the relevant parameters for the HL-LHC era and the FPF location. Above, we have also provided the relevant integrated luminosities. The multiple collisions that occur in each bunch crossing (pile-up) are accounted for in determining the flux of V . Throughout the study, we assume perfect detection efficiency for all the events that pass the selection criteria. The probability of passing such criteria depends on the geometrical acceptance of the detectors, energy and other kinematic cuts, as well as on the final state interactions inside the nucleus that we take into account in the case of the elastic scattering. We discuss the relevant cuts for different signatures below.
We will also include in our plots the expected sensitivities of the SND@LHC detector [98] to DM scattering in the U(1) B model, as determined in Ref. [99]. For the elastic DM scattering signature, this analysis assumed that backgrounds from muon-induced hadrons and photons can be rejected and that the number of neutrinoinduced events can also be suppressed to a negligible level. In the DIS regime, the analysis estimated that pure neutrino-induced backgrounds could be reduced to O(1000) events, and the sensitivity curves were taken to be N = 100 DM signal event contours.

IV. SIGNATURES
The hadrophilic models we are considering produce a diverse array of new physics signatures. These are shown in Table I, where we list which models are relevant for each signature, the dominant production and detection processes that determine the signal rates, the dependence of these rates on the model parameters, and the dominant SM backgrounds. As can be seen, the FPF experiments will be sensitive to direct signals generated by both the dark vector boson and DM, as well as to neutrino-induced signals. We now discuss them all in detail.

A. DM Deep Inelastic Scattering
We first consider DM DIS off nuclei, χN → χX. At large momentum transfer, DM DIS produces a significant hadronic recoil with multiple charged tracks. The main background is SM neutrino NC interactions. Due to the light mediator, DM scattering prefers lower momentum transfer than the neutrino background, which proceeds through Z-boson exchange. Our discussion of this signature closely follows that in Ref. [8].
The differential cross section for complex scalar DM DIS in the models of Sec. II is given by where 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, Q 2 = 2m p E χ xy is the squared momentum transfer, and f q is the quark parton distribution function. We use the nCTEQ15 parton distribution functions [100] for tungsten and argon and integrate Eq. (16) requiring Q 2 > 1 GeV 2 to obtain the expected numbers of DM DIS events in the FPF detectors. We also require the energy transferred to the hadronic system to be 1 GeV < E had < 15 GeV, where E had = yE χ , and the total transverse momentum of the recoiling hadrons to be 1 GeV < p T,had < 1.5 GeV, where p 2 T,had = Q 2 (1 − y). For the background, we calculate the expected numbers of neutrino NC scattering events satisfying the same cuts on Q 2 , E had , and p 2 T,had . Our cuts favor softer hadronic recoils, eliminating much of the neutrino NC background. Our projected sensitivities assume perfect detector efficiency and consider only statistical uncertainties. A previous study [8] of DM DIS at FLArE found that some experimentally motivated cuts did not have a large effect on the signal, but a full study remains to be performed.

B. DM-Nucleon Elastic Scattering
The light DM particles produced in the far-forward region at the LHC can also be discovered via their elastic scatterings with nucleons, which lead to single proton tracks visible in the detector. We treat this signature following Ref. [8], in which we have also studied the relevant neutrino-induced backgrounds. In particular, when presenting the sensitivity contours, we require the momentum of the outgoing proton to be within the range 300 MeV < p p < 1 GeV in FASERν2, and for FLArE we require p p < 1 GeV and the proton's kinetic energy to 8 Sec. IV A Sec. IV B Sec. IV C Sec. IV D Sec. IV E Models Rate scales as In the first three rows, the name of the signature, the subsection in which it is discussed, and the relevant new physics models are given. In the 4th and 5th rows, we show the Feynman diagrams for some example production and detection processes, respectively. The production processes shown are not necessarily the dominant ones. The 6th row shows the dependence of the signal rate on the model parameters, and the 7th row lists the dominant SM backgrounds.
satisfy E k,p > 20 MeV. We also reject events in which other visible tracks emerge from the vertex. After these cuts, we expect ∼ 100, 1000, and 300 background events during the entire HL-LHC run for FLArE-10, FLArE-100, and FASERν2, respectively. The elastic scattering cross section for the complex scalar DM interacting with the neutron or proton via the hadrophilic gauge boson is is the squared fourmomentum transfer in terms of the nucleon mass m N (N = n, p) and the outgoing nucleon energy E N , and E χ corresponds to the incident DM energy. The term proportional to A(Q 2 ), which contributes negligibly to the cross section at high energies, is given by with τ = Q 2 /(4m 2 p ). In contrast to the case of a vanilla dark photon mediator, the neutron and proton form factors are identical in this case and given bỹ where µ p = 2.793, µ n = −1.913, and G D (Q 2 ) = (1 + Q 2 /M 2 ) −2 , with M = 0.843 GeV. The differential elastic scattering cross section becomes form-factor suppressed at large momentum transfers, and the total elastic cross section is dominated by the contribution from Q 2 m 2 V . In the following, we include scatterings off both protons and neutrons. For protons, we include the efficiency factors ∼ (50 − 70)% related to the final-state interactions (FSI) of protons, as in Ref. [8]. For neutrons, we include similar efficiency factors in the range (15 − 30)%, which have been obtained as a function of the outgoing neutron momentum by studying neutrino interactions in GENIE [101,102]. In this case, the neutron re-scatterings inside the nucleus can lead to an outgoing proton with momentum within the aforementioned cuts and with no other detectable tracks. We find that scatterings of DM off neutrons can contribute up to 25% to the total elastic event rate.

C. Enhanced Neutrino Neutral Current Scattering
When a new mediator couples to neutrinos, NC scattering νN → νX receives an additional contribution from the mediator. The signature is identical to that for DM DIS. However, as NC scattering depends only on the couplings of the mediator to quarks and neutrinos, there is no dependence on m χ or Q χ , unlike the case of DM scattering. In particular, for the B − 3τ mediator, the total ν τ NC cross section becomes where Here g W is the SM weak coupling, g ν,L = 1 2 , and g q,L = 1 2 − 2 3 sin 2 θ W for up-type quarks and − 1 2 + 1 3 sin 2 θ W for down-type quarks. The second term in c L,R is the contribution from the new B − 3τ mediator with charges Q ν = −3 (3) for ν τ (ν τ ), and Q q = 1 3 for all quarks. The interference term is proportional to Q ν Q q , and so carries opposite signs for ν τ and ν τ NC scattering [103]. At the FPF where we expect almost equal fluxes of ν τ andν τ , this implies a small contribution from the interference term after cancellations. Nevertheless we use the complete expression above in our analysis.
For small m V , the BSM contribution to NC scattering prefers low recoil energy, similar to DM DIS and unlike the weak boson-mediated SM process, whose cross section grows with momentum transfer. We calculate the number of additional NC events expected at the FPF with Eq. (21), using the same parton distribution functions and minimum Q 2 cut as in Sec. IV A. Because of the small relative flux of tau neutrinos compared to muon and electron neutrinos, the impact of the light mediator on the total NC cross section must be significant to provide a sizable effect relative to the SM NC background.
In testing whether an excess of NC events is observable, we consider only statistical uncertainties and neglect systematic uncertainties. For simplicity, we also assume perfect detection efficiency for NC interactions; the inclusion of realistic detection efficiencies [104] would not substantially change the positions of the limits from excess NC events in Figs. 2 and 3, relative to the other signatures that we consider. We note that the main systematic uncertainty in the NC cross section measurement, the neutrino flux, can be constrained by measurements of charged current (CC) interactions. We find a statistically significant effect from the BSM contribution to NC scattering when the coupling to mass ratio of the new interaction is comparable to that of the weak interaction,

D. Excess of Tau Neutrino Flux
In the case of the gauged B − 3τ scenario, the hadrophilic mediator decays into tau neutrinos with a sizable branching fraction. 5 As discussed in Ref. [51], this opens another opportunity to probe this model via their contribution to the LHC tau neutrino flux. In the SM, tau neutrinos are mainly produced via D s → ντ and subsequent τ decays, which occurs in roughly one in 10 5 collisions at the LHC. This means that even rare BSM processes could lead to sizable contributions to the tau neutrino flux. The relevant detection channel in this case is via ν τ CC scatterings off nuclei. The displaced decays of the outgoing boosted tau lepton must then be identified in the detector, requiring excellent spatial resolution.
An important issue that arises when searching for signs of new physics is the large uncertainty on the normalization of the SM tau neutrino flux [96,105]. Although future efforts are expected to reduce these uncertainties, we will follow a different approach. In contrast to tau neutrinos from charm and tau decays, which have a broader angular spread, tau neutrinos from light mediator decays are more centered around the beam collision axis. In this study, we use this feature and perform a shape analysis of the ν τ angular distribution, which does not rely on knowledge of the neutrino flux normalization. We focus on the FLArE-10 design, whose 1 m × 1 m cross sectional area is sufficiently large to capture this effect. More precisely, we define five concentric rectangular bins centered around the beam collision axis and corresponding to the distance between d and d + 10 cm away from it, where d = 0, 10, 20, 30, 40 cm. In practice, the most important contribution to the BSM-induced excess of ν τ s is from the two most central bins, i.e., at distances up to d 20 cm away from the beam collision axis.

E. Visible Decays of the Dark Vector Boson
In the following, we study the decay signature using the FORESEE package [41] with the lifetimes modeled with DarkCast [45] and the spectrum of light hadrons obtained from EPOS-LHC [106]. We assume 100% detection efficiency for all visible final states. We present the results for both FASER and FASER2. In the analysis, we require the total energy of the visible products of the vector boson decays to be at least 100 GeV. This cut has a minor impact on the BSM signal events, but suppresses possible SM backgrounds to a negligible level [24,97]. Visibly decaying dark vector bosons could also appear in secondary production processes due to DM scatterings occurring right in front of or even inside the detector [107]. We neglect the impact of such processes below, as we do not expect them to improve the sensitivity reach of the FPF detectors in the models under study.

V. RESULTS
In Fig. 2, we present the results of our analysis for both the U(1) B and U(1) B−3τ models in the (m V , g V ) plane. In the plots, we fix the DM coupling to α χ = 0.01 and 0.5 in the upper and lower panels, respectively, and we keep a constant mass ratio between the dark sector particles, m V = 3 m χ . In dark gray, we show the existing constraints, as discussed in Sec. II D, while the black solid (dashed) lines correspond to the relic density targets for the complex scalar (Majorana) DM. We stress that, although the anomaly bounds, shown in light gray in the left panels for the U(1) B case, can be avoided in modified versions of this simplified scenario, this often leads to further constraints due to additional couplings of the dark vector bosons that are introduced in the model to make it anomaly-free. An example is shown in the right panels for the anomaly-free U(1) B−3τ model, where the NSI constraints cover a good portion of the parameter space shown in the plot.
In Fig. 2, we also present the expected 90% CL exclusion bounds in searches for DM scatterings off nuclei in the elastic (dark red) and DIS (light red) channels for FLArE-10 (solid), FLArE-100 (dash-dotted), and FASERν2 (dotted). As is clear from the plot, the elastic scattering probe is stronger for light DM and mediator masses below 1 GeV, which favor interactions with low momentum exchange. For m V 1 GeV, the elastic scattering rate is suppressed by the form factor and the cut on the outgoing proton momentum p p 1 GeV. In this higher mass range, the search based on DIS processes provides the best reach. For comparison, we also show the expected reach of the SND@LHC detector [99] with the assumptions noted in Sec. III.
For the U(1) B model with fixed α χ = 0.01 shown in the upper left panel of Fig. 2, we expect that future searches at the FPF will cover almost the entire remaining allowed region in the parameter space above the Majorana and complex scalar relic target lines, in which DM is not thermally overproduced in the early Universe. This corresponds to vector boson masses between 1 and 3 GeV. For the simple complex scalar DM model, additional stringent bounds for m χ 200 MeV can arise from past DM DD searches, which are indicated in the plots by the very light gray shaded regions and cover the region within the sensitivity of FLArE and FASERν2. However, these limits can be evaded in the inelastic scalar DM case and are not relevant for Majorana DM. For lower masses, (a few) MeV m V 1 GeV, the expected FLArE and FASERν2 bounds extend beyond current constraints from the CCM, MiniBooNE, and NA62 experiments. Here, the searches at the FPF would probe regions in the parameter space that are otherwise partially excluded only by anomaly-induced rare K and Z decays.
Next, we consider the U(1) B−3τ model with fixed α χ = 0.01 shown in the upper right panel of Fig. 2. Since the model is free of gauge anomalies, the stringent con-straints from rare Z and meson decays present in the U(1) B model are absent in this case. On the other hand, the additional bounds from neutrino NSI cover much of the model parameter space. Nevertheless, we observe that the FPF detectors can still explore a portion of the currently allowed parameter space, especially in the ω and φ resonance regions, m V ∼ m ω , m φ , and the corresponding part of the relic target line for complex scalar DM. In this model, additional sensitivity arises from dark vector boson-mediated scattering of tau neutrinos in the DIS regime; see Sec. IV C. The relevant expected bounds, which are indicated by the light purple lines in the plots, impact parameter regions that are already excluded by past searches. We note that the actual exclusion bound in the DIS channel should be derived using the combined excess signal rates for both DM and BSM neutrino scatterings over the expected SM backgrounds. Instead, in the plot, we have presented the expected bounds for each separately to allow for independent discussion of the impact of different new physics effects.
For larger values of α χ , the relic target lines shift downwards relative to the FPF sensitivity contours from DM scattering. This is dictated by the different parametric dependence of the annihilation cross section and the number of DM scattering events in the FPF experiments on the coupling constants, σv ∼ g 2 V α χ and N ev ∼ g 4 V α χ , respectively. As a result, in the lower panel of Fig. 2 obtained for α χ = 0.5, we observe that both FLArE and FASERν2 will only partially cover the thermal target lines for the U(1) B model. Instead, in the U(1) B−3τ case, they will typically probe regions in the parameter space predicting subdominant fractions of thermally-produced χ DM.
Thus far we have considered scenarios in which the vector boson mediator couples much more strongly to DM than to SM particles, Q χ 1. In Fig. 3 we consider the different scenario in which the vector boson mediator couples with comparable strength to complex scalar DM and SM particles, with Q χ fixed according to Eq. (9). As can be seen, for both the U (1) B and U (1) B−3τ models, FLArE-10 can cover the entire relic target line in a wide vector boson mass range between 1 MeV and 10 GeV. As in the previous scenarios depicted in Fig. 2, significant portions of these regions are already constrained by either anomaly-induced or NSI bounds, as well as by the other past searches indicated in the plots. However, we emphasize that for the case of inelastic scalar DM in the U(1) B model, to which DD constraints do not apply, FLArE-10 will be able to test an interesting open region of parameter space for vector boson mass of order several GeV that is consistent with the observed DM abundance. DM is thermally overproduced below these contours. The light (dark) red lines correspond to 90% CL exclusion bounds from DM DIS (elastic) scatterings off nuclei for FLArE-10, FLArE-100, and FASERν2, as indicated. The dotted brown contours are the sensitivity contours for SND@LHC [99]. In the right panels, the light purple contours are the projected sensitivity contours from probing the V -induced BSM NC interactions of tau neutrinos. In both panels, the dark gray shaded regions are excluded by current bounds. The light gray shaded regions in the left (right) panels correspond to the anomaly-induced K and Z decays (NSI bounds). The very light gray shaded regions are constraints from DM DD; these do not apply to Majorana and inelastic scalar DM (see Sec. II D).
between several hundred MeV and a GeV and coupling 10 −8 g V 10 −5 , such displaced decays into visible final states, primarily light hadrons, can be detected at both FASER and FASER2. The dominant branching fraction in this case is into three pions, π 0 π + π − , which leads to a striking signature consisting of a photon pair accompanied by two oppositely charged tracks.
We present the relevant expected 90% CL exclusion bounds on LLP decays for FASER (dark blue) and FASER2 (light blue) in the plots. These correspond to the region of parameter space with m V ∼ m ω or m φ . Here, both the dark vector boson production via proton bremsstrahlung and its decay branching fractions into light hadrons are enhanced. The expected exclusions shown in the plot are bounded from below by the production rate of the dark vector bosons being too low, and from above by the V lifetime being too small for the boson to decay in the detectors. In the U(1) B model, further sensitivity at FPF experiments can be obtained for m V 10 MeV due to loop-induced dark vector boson decay into an e + e − pair. This scenario is, however, already constrained by the past beam-dump search in NuCal and by the anomaly-induced bounds.
Last but not least, in the U(1) B−3τ model, further constraints arise due to the dominant dark vector boson decays into tau neutrinos. These can generate an excess flux of ν τ s over the expected SM production rates, which can be detected via their CC scatterings in the detector, Same as Fig. 2, but for only the FLArE-10 detector, complex scalar DM, and fixed charges Qχ = 1 (left) and Qχ = −Qτ = 3 (right), resulting in a floating αχ. Additional expected exclusion bounds from probing displaced V decays to SM final states in FASER (FASER2) are shown with dark (light) blue lines. In the right panel, the green contour is the sensitivity contour from probing excess CC scatterings of ντ .
as described in Sec. IV D. The corresponding expected sensitivity is indicated by the green contour in the right panel of Fig. 3. For m V 2 GeV, this sensitivity is greater than from the DM and BSM neutrino searches. In particular, it allows one to constrain the currently allowed region of the parameter space of the model close to the ω-and φ-resonance regions. In this case, the increased flux of ν τ s could also further contribute to the aforementioned NC DIS signal rate due to BSM tau neutrino interactions. To isolate the impact of various new physics effects, we do not take this into account when presenting relevant expected bounds, which should thus be considered conservative. We stress that the dominant expected bound in the corresponding region of the parameter space is, in any case, due to excess CC ν τ scatterings.

VI. CONCLUSIONS
While the ability of the LHC to search for TeV-scale DM is well known, recently proposed dedicated experiments at high rapidity can significantly enhance the potential of the LHC to probe light DM. Beyond the mini-mal portal extensions of the SM that allow for light DM and an associated mediator, new gauge groups represent a well-motivated class of possible dark sector models. In this paper, we have explored the use of the FPF to study such U(1) theories leading to hadrophilic dark sectors. In particular, these remain beyond the reach of experiments focusing on BSM electron couplings, while they can more straightforwardly be studied at the LHC.
The suite of FPF experiments provides a comprehensive set of tests of these theories in different regions of parameter space. DM produced in pp collisions can scatter in the FASERν2 and FLArE detectors through the new light vector boson; we have considered both elastic and deep inelastic scattering. Furthermore, if the mediator has a significant decay branching ratio to SM states, the FASER2 LLP detector can search for the visible decay products. In fact, already at Run 3, the FASER detector will begin to test hadrophilic U(1) theories at couplings substantially lower than existing bounds. These hadrophilic models therefore motivate near-term searches at FASER for new LLP signatures V → π 0 γ, π + π − π 0 , K + K − , K S K L , which are not motivated by dark photon models for 13 FASER in Run 3. Finally, if neutrinos are also charged under the new gauge symmetry, additional signatures are possible in the scattering detectors. We have demonstrated that with a symmetry under which tau neutrinos are charged, the ν τ flux and NC cross section are both enhanced, leading to potential deviations in ν τ CC and NC scattering rates.
These results for U(1) B and U(1) B−3τ models should be considered as illustrative of the complementarity of forward LHC experiments in searching for light dark sectors, particularly between LLP and scattering detectors. In both of these theories, the FPF can test broad regions in the coupling-gauge boson mass plane, including significant expanses over which the observed DM relic density could be obtained through standard thermal freezeout. For our benchmark scenario with scalar DM, m V = 3m χ and low values of the dark charge Q χ given by Eq. (9), FPF searches can probe well below the thermal relic target lines in each model for nearly all gauge boson masses between 1 MeV and 10 GeV. Throughout our results, the strongest searches tend to be those based on DM elastic scattering and DIS, with distinct additional reach possible from LLP searches when the mediator can decay to SM final states. For the B − 3τ model, searches for an increased ν τ flux would also test new space.
The models that we have studied face strong indirect constraints, notably from rare invisible decays and neutrino oscillations, but we emphasize that FPF searches can test couplings that are smaller than these formidable existing bounds. In addition, in the GeV mass range, these searches provide constraints that are complementary to those from spin-independent DD, the latter only being applicable in the case of elastic scalar DM.
Though we have chosen to focus on two possible gauge groups with a handful of coupling and mass assumptions, the general interplay between the DM scattering, LLP and neutrino searches is likely to persist for other theories and parameter choices. To determine the gain provided by the FPF in a particular theory, the reach of these searches must be compared against those from other bounds. As we have seen, U(1) theories that are not anomaly-free typically face rare meson decay constraints, while those with nonzero lepton charges can encounter NSI bounds. For models with couplings to 1st and 2nd generation leptons, additional limits from beam dump and neutrino experiments would likely need to be considered as well.
Forward LHC detectors offer a distinct perspective on light hidden sectors, allowing for searches for light DM and its associated mediators in an otherwise inaccessible kinematic regime. The results here underscore the utility of different types of forward detectors, as could be provided at the FPF. The multi-pronged approach to uncovering physics beyond the SM that is enabled by such a facility, along with other uses such as measurements of SM neutrino interactions and tests of QCD, bolsters the physics case for the FPF.