FLArE up dark sectors with EM form factors at the LHC Forward Physics Facility

Despite being mostly secluded, dark sector particles may feebly interact with photons via a small mass-dimension 4 millicharge, a mass-dimension 5 magnetic and electric dipole moment, or a mass-dimension 6 anapole moment and charge radius. If sufficiently light, the LHC may produce an intense and collimated beam of these particles in the far forward direction. We study the prospects of searching for such dark sector particles with electromagnetic form factors via their electron scattering signature in the Forward Liquid Argon Experiment (FLArE) detector at the Forward Physics Facility (FPF). We find that FLArE can provide new probes of sub-GeV dark particles with dipole moments and strong sensitivities for millicharged particles in the 100 MeV to 100 GeV region. This complements other search strategies using scintillation signatures or dark matter direct detection and allows for probing strongly interacting dark matter motivated by the EDGES anomaly. Along with the FORMOSA detector, this leads to a very diverse and leading experimental program in the search for millicharged particles in the FPF.


I. INTRODUCTION
Shedding light on beyond the Standard Model (BSM) constituents of nature remains of primary importance for our understanding of elementary interactions. Of particular interest are possible couplings of the new species to the Standard Model (SM) photons that could leave an impact on the entire cosmological history, as well as could nowadays be probed experimentally. In particular, while dark matter (DM) species are generically expected to have suppressed such interactions, even very weak couplings of this type could lead to striking observational consequences. Therefore, pushing the boundaries of our understanding of BSM couplings to photons helps us understand how the dark sector of the Universe works.
Such couplings are generally expected to occur in many BSM scenarios, either directly or at a loop-induced level. At low energies, the resulting couplings can be written in terms of effective field theory (EFT) operators, which generate electromagnetic (EM) form factors of the dark species. This can lead to both electrically charged and neutral new particles with suppressed interactions with the SM photons that can be within reach of current and future searches.
The far-forward region of the LHC presents new exciting opportunities to conduct novel searches of such BSM species. At this unique location, even distant and small detectors can search for new light unstable particles [71][72][73][74], as well as study high-energy neutrinos [75][76][77][78]. This observation lead to the approval of FASER [79][80][81][82][83][84], FASERν [85,86], and SND@LHC [87,88] detectors to take data during LHC Run 3, as well as to the proposal of a dedicated Forward Physics Facility (FPF) [89][90][91] for the High Luminosity LHC (HL-LHC) era. The FPF would host several experiments, including the proposed FORMOSA detector [12] in which mCPs could be searched for via a low ionization signal. Compared to other mCP studies, this search benefits from the TeVscale center-of-mass (CM) energy of pp collisions at the LHC such that mCPs up to the mass of order few tens of GeV can be probed with O(10 −3 ) EM charge. The focus on the far-forward direction allows for maximizing the flux of mCPs such that a relatively compact detector can provide strong bounds on this scenario in the future.
In this paper, we discuss alternative detection signatures based on scatterings of high-energy dark particles with EM form factors producing detectable soft electron recoils in the FPF detectors. To this end, we present our results, for the recently proposed Forward Liquid Argon Experiment (FLArE) experiment [92] which is envisioned to employ liquid argon time projection chamber (LArTPC) technology to detect visible signals from such scatterings, and to study neutrino interactions [90,91]. Given the excellent capabilities of LArTPC detectors to study soft electron-induced signals at the energy level of above O(10 MeV), and to reconstruct such events, FLArE will be well-suited to study low-energy scattering signatures. We also stress, however, that similar sensitivity could be expected for other FPF scattering experiments, i.e., Advanced SND@LHC and FASERν2 detectors, cf. Refs [90,91] for further discussion, provided their final design allows for sensitive searches for low-energy electron recoils.
As shown below, the combined research agenda of the FPF will lead to a variety of experimental approaches and world-leading detection prospects in the search for mCPs in a wide range of their masses, with potentially important connections to the aforementioned astrophysical observations. We also present, for the first time, expected FPF bounds for the higher-dimensional operators. In particular, we show the projected bounds for electric and magnetic dipole moment that can improve current constraints for dark states in the sub-GeV mass range.
This paper is organized as follows. In Sec. II, we describe the considered particle model. In Sec. III, we discuss dominant dark state production channels, followed by detection strategy and background estimation in Sec. IV. In Sec. V, we demonstrate the projected sensitivity of FLArE and compare it with existing constraints. Finally, we draw a conclusion in Sec. VI.

II. DARK STATES WITH ELECTROMAGNETIC FORM FACTORS
We consider a scenario in which the SM photon is the new physics mediator which couples to dark sector particles χ carrying EM form factors. To have the most affluent set of EM form factors, in this work we take χ to be a Dirac particle. Starting from the lowest dimensions and till mass-dimension 6 in the EFT framework, 1 allowed effective operators are mass-dimension 4 electric monopole, mass-dimension 5 magnetic/electric dipole moment (MDM/EDM) and mass-dimension 6 anapole moment (AM) and charge radius (CR). The Lagrangian describing interactions between χ and the SM photon can be expressed as [57], where F µν is the field strength of the SM photon, is the fraction of χ's charge relative to the elementary charge e, µ χ and d χ are coefficients of MDM and EDM with dimension [M ] −1 , a χ and b χ are coefficients of AM and CR with dimension [M ] −2 , and σ µν ≡ i[γ µ , γ ν ]/2. Below the electroweak scale, the effective couplings to photons in Eq. (1) can be induced from more general hypercharge couplings. At mass-dimension 4, the tiny millicharge can generally arise from a small hypercharge e carried by χ before the electroweak symmetry breaking. The interaction Lagrangian between millicharged χ and the SM U(1) Y gauge field B µ can be expressed as eχγ µ χB µ . This term can be added to explicitly break the charge quantization, or can be more theoretically motivated and source from a kinetic mixing between a dark photon and the SM photon [5]. After electroweak symmetry breaking, we then observe that = cos θ W with θ W being the weak mixing angle. Similarly, the coupling between χ and the Z boson will be induced. We discuss its impact on the sensitivity reach of the FPF experiments for the mCP model in Sec. V.
In the case of higher-dimensional operators, similar issues arise from replacing F µν with the hypercharge field strength tensor B µν in the operators in Eq. (1). Notably, this has important consequences for both the missing energy searches for χs at the LHC and its DM phenomenology for m χ 100 GeV, as in both cases one probes regions of the parameter space beyond the EFT validity regime for the couplings to F µν [59]. Instead, at lower energy scales characteristic for light hadron decays, the expected sensitivity is driven by the photon couplings. We then find it sufficient to focus on this regime below, while we also comment on consequences of employing the hypercharge couplings.
Once the hypercharge coupling is taken into account, the interaction is a priori UV-completed for massdimension 4. In contrast, for mass-dimension 5 and 6 operators, EFT is not valid anymore when we start to probe the full particle spectrum of the underlying theory. The effective theory considered in this work can be UV-completed by considering, e.g., compositeness of χ [93][94][95] or loop contribution via new electric-charged states at UV-scale [96,97]. Assuming the UV-scale is not within the reach of FPF, we adopt the effective Lagrangian Eq. (1) in the following calculations.

III. PRODUCTION CHANNELS AT THE FORWARD PHYSICS FACILITY
The LHC is the most energetic particle collider build thus far and one of its primary objects for its remaining operations is the search for new heavy particles at the TeV scale. These particles are expected to decay roughly isotropically, and large experiments have been build around the LHC interaction points to detect their decay products. As first pointed out in Ref. [71], the LHC could also be used to probe a different class of new particles which are light and weakly coupled particles. This idea utilizes that the LHC produces an enormous number of hadrons, which are mainly produced in the direction of the colliding beams and typically escape undetected through the beam pipe. If kinematically allowed, these hadrons could undergo rare decays into new light weakly coupled states and produce an energetic and strongly collimated beam of these particles in the forward direction. In this study, we perform a dedicated Monte Carlo simulation using FORESEE [98] to estimate the flux of dark sector particles χ produced at the LHC.
If sufficiently light, χ can be produced through both three-body decays of pseudoscalar mesons P → γχχ and two-body decays of vector mesons V → χχ. In the following, we consider the pseudoscalar mesons P = π 0 , η, η as well as the vector mesons V = ω, ρ, φ, J/Ψ, Ψ(2S), Υ(nS) as parent particles and use particle spectra provided by FORESEE [98]. Here, the spectra of light mesons were obtained using the EPOS-LHC [99] event generator as im-plemented in the CRMC [100] interface. For the heavy Charmonium and Bottomonium states, we use the spectra presented in Ref. [12] which were obtained using Pythia 8 [101,102] and tuned to LHCb data [103][104][105].
Heavier χ, for which the above-mentioned meson decays are kinematically forbidden, are primarily produced via the Drell-Yan process qq → χχ. We simulate this using MadGraph 5 [106] with the model file adapted from Ref. [65]. To ensure that the parton distribution functions are well defined, we only consider this partonic scattering production mode for masses m χ > 1 GeV. For smaller masses, and as a result also smaller momentum transfers, the partonic picture loses its validity. In this case, other processes such as production via coherent Bremsstrahlung off the proton beam could play a role [107].
In Fig. 1 we present the production rates of the dark sector state χ at the LHC with 14 TeV center-of-mass energy for the different production channels mentioned above. The left panels correspond to the total production rate, while the right panels show the production rate after requiring the particle to be within 1 mrad around the beam collision axis, which roughly corresponds to the angular size of FLArE.
For mCPs, the production rate is roughly energyindependent. Thus χ-production is dominated by decay of lighter mesons due to their larger fluxes. However, the hierarchy is different for higher-dimensional operators for which production rate is energy-dependent. For example in the case of two-body vector meson decays, their decay branching fractions scale as BR V →χχ ∼ 1 for mCPs, BR V →χχ ∼ M 2 for EDM/MDM and BR V →χχ ∼ M 4 for CR/AM, where M is the vector meson mass (see Appendix A for more details). As a result, the importance of heavier meson decays or the Drell-Yan (DY) production increases for dark sector states with higher dimensional operators. After imposing the angular cut, we find that for mass-dimension 5 operators (MDM and EDM) production from J/ψ decay can be comparable to lighter meson decay, while for mass-dimension 6 operators (AM and CR) production from J/ψ decay and the DY processes remain dominant throughout the considered m χrange. Notably, while the DY production can probe large center-of-mass energies in the hard collision, for which mass-dimension 5 and 6 operators with only an F µν coupling could suffer from the aforementioned lack of the EFT validity, we will see below that this production process does not determine the FLArE sensitivity reach in the currently allowed regions of the parameter space of these models.
In the case of the mCP model, we also recognize in the upper panels of Fig. 1 the impact of the Z boson coupling induced by small hypercharge. This can be seen as a characteristic peak of the total DY production rate for m χ m Z /2, where the effective mCP production via the Z resonance occurs. We expect similar peaks to be present for other models if the hypercharge coupling was taken into account (not shown in the plots). This, how- FIG. 1. Production rates of dark states with millicharge (top panels) and EM form factors with MDM (center) and CR (bottom) within the whole forward hemisphere (left panels) and within 1 mrad of the Line of Sight (LoS) (shown on the right). Due to the nature of effective operators, production from heavier mesons and the Drell-Yan (DY) process becomes increasingly important for higher mass-dimension, in contrast to the mCP case where production from π 0 decay dominates; see main texts for further details. Note that production rates of EDM (AM) case only differs with that of MDM (CR) at the kinematic endpoints, thus we do not show them here. ever, does not affect lighter χs for which the dominant production is through hadron decays.

IV. SIGNATURES AT FLARE
Stable dark sector particles which are produced in the far-forward region of the LHC can be studied via their scatterings in the detectors in the FPF. In our discus-sion, we focus on electron-recoil signals induced by new physics species. These can be more efficiently disentangled from neutrino-induced backgrounds than signatures based on scatterings off nuclei [108,109], especially in searches targeting low energy depositions in the detector. In the following, we discuss several strategies to maximize the sensitivity reach and we comment on expected backgrounds.

A. Signal
Scattering a-la DM signal -The first signature that we consider resembles the single electron-recoil signature, which has previously been considered for a light DM search in FLArE [92]. In the relativistic regime, we find for the scattering cross section of the process χe → χe that dσ/dE R ∝ E x R where E R is the electron recoil energy and x = −2, −1, 0 for mCP, MDM/EDM and AM/CR, respectively; see Eq. (B2). Therefore we can infer that for lower mass-dimension operators, the event rate is enhanced for softer recoils, which is especially important for mCP searches [12]; see also the following discussion.
In our analysis below, we assume the following energy thresholds for the recoiled electron where the lower cut of 30 MeV corresponds to the assumed LArTPC detector capabilities to study soft electron tracks, cf. Refs [110]. In order to illustrate the impact of this cut on our analysis, we also present below the expected FLArE sensitivity reach in the search for mCPs assuming an increased lower energy threshold which is set at 300 MeV. The upper energy threshold of 1 GeV allows for suppressing neutrino-induced backgrounds while maintaining the signal rate, as discussed below.
Double-hit with softer recoils -Since the differential scattering cross section of mCPs is more IR-biased, further improvement in the sensitivity reach can be expected for an even lower energy threshold in LArTPC detector [37]. While lowering it too much could lead to additional backgrounds, these can be circumvented by focusing on double-hit events, in which the signature consists of two coincident and collinear χ scatterings off electrons [40]. We, therefore, present the expected sensitivity of FLArE to such a search for mCPs assuming the following cuts on both electron recoil energies where the lower energy threshold to detect each of the hits in argon is chosen following Ref. [110]. The mCPinduced signal, in this case, would consist of two simultaneous hits which, given the large boost factors and very small deflection angles of χs, define the line pointing towards the direction of the ATLAS IP.

B. Background
The experimental signature of our interest in FLArE consists of a single scattered electron or two coincident such electron recoil signals. The signatures of this kind can also be mimicked by other types of effects. This leads to possible backgrounds that we briefly discuss below. We first focus on backgrounds induced by the two types of SM species which can reach the FPF experiments while being produced at the ATLAS IP, namely neutrinos and muons. We then comment on possible other sources of backgrounds that remain more difficult to estimate without detailed detector simulations.
Neutrino-induced backgrounds -Signatures of new physics particles scattering in the detector often resemble similar processes characteristic to the SM neutrinos. In particular, the signal consisting of single scattered electrons in FLArE has already been studied in the context of light DM searches in Ref. [92]. It has been shown that for low visible energy depositions in the detector, 30 MeV E R,single,loose 20 GeV, and thanks to the use of additional angular cuts, the expected neutrinoinduced background rate can be as low as O(20) events during the HL-LHC era in the 10-tonne FLArE detector placed in the FPF along the beam collision axis. The most significant background contribution is associated with the neutrino-electron scatterings, while quasielastic and resonant nuclear scatterings of ν e also contribute non-negligibly. Instead, neutrino-induced deep inelastic scattering events can typically be rejected due to the presence of additional visible tracks in the detector.
Notably, the number of neutrino-induced events expected in the FPF is highly suppressed in this low-energy regime. This is primarily dictated by the large energy of incident neutrinos produced in the far-forward region of the LHC, typically of order E ν ∼ 200 − 300 GeV, see Refs [90,91,[111][112][113][114] for further discussion. Instead, low-energy neutrinos are produced more isotropically at the LHC. As a result, the remaining low-energy background from ν − e scatterings in FLArE is primarily associated with interactions of high-energy neutrinos with E ν > 100 GeV that can, occasionally, generate soft electron recoils.
We employ this fact in the analysis below by even further reducing the electron recoil energies required for our signal events in the single scattered electron signature, cf. Eq. (2). In this case, we find less than O(1) expected neutrino-induced background events in FLArE, while the number of coincident double scattering events is significantly lower. In our estimates, we use far-forward neutrino fluxes and spectra presented in Ref. [108] obtained using the fast neutrino flux simulation introduced in Ref. [112] and we use GENIE [115,116] to study neutrino interactions.
Neutrino-induced backgrounds can also appear due to ν scatterings in the rock and other material in front of FLArE. Such interactions can induce secondary particles and EM showers entering the detector from outside. Charged such species can be vetoed, similarly to the vetoing capabilities envisioned for FASER/FASER 2 detectors [79,81]. Instead, rare events with secondary neutral hadrons entering the detector's fiducial volume with no charged counterpart reaching FLArE would be more difficult to veto. Such events could mimic single scattered electron signatures by inducing additional photons in the detector, if these are energetic enough and if the remain-ing neutral hadrons remain undetected. We leave a detailed analysis of this background contribution for the future detector simulation. Since they could correspond to only a tiny fraction of neutrino interactions in front of the detector, we assume these background remains suppressed in the current analysis.
Muon-induced backgrounds -Instead, high-energy muons passing through the FPF detectors could more straightforwardly generate signals mimicking scattered electrons. This is primarily due to muon-induced photons produced in the detector [92]. In order to suppress the rate of muons traveling through the FPF, the installation of a dedicated sweeper magnet has been proposed in the LHC tunnel which could deflect forward-going muons before they reach the experimental facility [90,91]. In addition, such events are expected to be vetoed by detecting a charged muon entering FLArE. For this purpose, in order to improve event timing and triggering beyond the limitations of the drift time in LArTPC detectors, it is envisioned to equip FLArE with an additional light readout system [91]. Muons could also generate secondary charged and neutral hadrons due to their interactions both inside and outside of the detector, cf. discussion in Ref. [85] for the FASERν experiment. In the following, we assume that all such events are vetoed or rejected as not signal-like, as discussed above for neutrino-induced events.
Other background sources -Further backgrounds corresponding to empty detector time frames with no neutrino scatterings in the detector can appear, e.g., due to ambient gamma-ray activity, intrinsic radioactivity, or electronic noise, cf. Ref. [117] for discussion for the ArgoNeuT detector. In particular, these have been identified as important limitations of the mCP search in LAr detectors assuming order MeV detection thresholds [40]. Instead, as discussed above for the single scattered electron signature, we employ the increased threshold of order 30 MeV or so, which leads to proper EM showers detectable in LArTPCs, cf. e.g., Ref. [110] for further discussion. In the case of the double-hit signature, a significant background reduction can be achieved even for lower energy thresholds ∼ 5 MeV thanks to the requirement that both hits are aligned [40].

V. PROJECTED SENSITIVITY OF FLARE
We study the sensitivity reach for the FLArE detector placed in the FPF on the beam collision axis at a distance L = 620 m away from the ATLAS IP. We assume the relevant detector geometry (fiducial volume) where ∆ is the detector length along the beam collision axis and S T is its transverse area. In our study, we primarily focus on the 10-tonne FLArE detector, while for higher dimensional operators we also show how the expected reach could be improved for a larger 100-tonne experiment FLArE-100 with ∆ = 30 m and S T = (1.6 m × 1.6 m), cf. Ref. [92] for a similar comparison for the DM search. At this stage, we assume 100% signal detection efficiency for the dark states and study the electron scattering signatures discussed above.
In the plots below, we present the result corresponding to N ev = 3 expected new physics events in the detector. Assuming that backgrounds can be suppressed to a negligible level, this would correspond to sensitivity reach of FLArE in the proposed searches.

A. Millicharged particles and dark matter
We show the sensitivity reach of FLArE in the search for mCPs scattering off electrons in the detector in the left panel of Fig. 2, together with previous constraints shown as gray-shaded regions and projections of other future experiments. In the plot, the red solid line corresponds to a single-scattered electron and the recoil energy cuts of 30 MeV E R,single 1 GeV. As can be seen, at the low m χ end FLArE is able to compete with previous intensity frontier experiments such as LSND [37] and SLAC mQ [33], while for the mass range between 100 MeV and O(100 GeV) it can go beyond the past bounds from ArgoNeuT [41], BEBC [42], CMS [24,25], LEP [6], and the milliQan prototype [30].
When compared with future projections, although the FORMOSA detector [12] to operate in the FPF is the leading proposal based on its projection, the FLArE search based on electron scatterings could provide an independent measurement method for a wide range of masses for mCP with charges O(10 −3 -10 −1 ). It can also have better sensitivity in the mass range above ∼ 4 GeV than proposed FerMINI [38] and SUBMET [118,119] detectors. In the high-mass regime, FLArE is complementary to the sensitivity of milliQan [28][29][30][31], although it has a weaker sensitivity for m χ (40-50) GeV.
Importantly, while the other aforementioned proposed experiments will be targeting soft scintillation signals expected from mCPs traversing the detectors, the considered FLArE search is focused on larger energy deposits which can also be occasionally expected from energetic mCPs produced at the LHC. In order to highlight this even better, we also show in the plot with the red dashed curve the expected FLArE sensitivity assuming an increased lower recoil energy threshold of 300 MeV, cf. Eq. (2). In this case, the sensitivity reach becomes weaker, as expected for mCPs favoring soft electron recoils in the detector, cf. Eq. (B2) for the scattering cross section. However, it can still constrain important and currently allowed regions in the mCP parameter space.
Switching to multiple electron recoils, we find that the double-hit signal does not improve the sensitivity as the assumed detection threshold of 5 MeV is not low enough to identify smaller energy depositions characteristic for mCPs. This is in contrast to the scintillator-based de- 3 GeV) and the one with a lower energy threshold (ER 1 MeV) but requiring a double-hit event. The right panel shows the constraints for mCP dark matter, assuming 0.4% of the dark matter being mCP. One green band (0.4%) shows the parameter of mCP as dark matter to explain the EDGES anomaly, the other shows a model of 0.01% dark matter being mCP, while the mCP interact with the rest of the cold DM to achieve additional cooling; see main texts for details.
tector of FORMOSA, which is able to detect even O(eV) energy depositions. To similarly utilize the enhancement of event rate at low E R at FLArE, one needs a strategy that allows one to dial down the detection threshold, for example, a faint track signature. We leave this for future dedicated studies.
In the right panel of Fig. 2, we present a similar FLArE sensitivity due to single-and double-hit signals but assuming mCPs to be constituents of DM in the Universe. In this case, interesting effects on early-Universe cosmology and additional constraints from direct direction searches apply [9,120]. We present the latter with light gray-shaded regions in the plot assuming the mCP abundance to be 0.4 % of the total DM relic density to avoid strong cosmological constraints [21,121,122]. The parameter space we consider corresponds to the "strongly interacting dark matter" region [9][10][11][12], in the sense that the dark matter flux is attenuated through scattering in the atmosphere and crust. The mCP therefore loses its energy before getting to the terrestrial direct-detection experiments. Labeled as "Direct Detection" curve in the right panel of Fig. 2, we show such a curve of critical cross-section. Still, a dedicated balloon-based experiment dedicated to search for such strongly interacting dark matter [9] (labeled as RRS), and a rocket-based experiment, X-ray Quantum Calorimeter (XQC) [123] provide constraints above this curve as they can detect DM particles above the atmosphere (RRS and XQC data are recast as bounds in [10]). FLArE will provide a strong probe for the region above this direct detection curve, independent of the assumptions of DM abundance as mCP.
In addition, we also highlight slices of the parameter space which were used to explain the aforementioned EDGES anomaly [13,18,19,23,34] in Fig. 2. We show two interesting bands of parameter space. One light green band corresponds to a minimal scenario which mCP makes up 0.4% of dark matter, without additional interactions to the rest of the cold DM. We also plot two scenarios in which 0.01% (10 −4 ) of dark matter is mCP, and the mCP DM is coupled to the rest of the cold DM through a light mediator to achieve additional cooling of the gas to explain the anomalous absorption spectrum observed by the EDGES collaboration. FLArE, with both single-and double-hit considerations, can study this parameter space close to the EDGES predictions.

B. Dark states with electromagnetic form factors
For higher-dimensional EM form factor interactions, we focus on a single electron-recoil signal, since the scattering cross section does not favor low-energy deposition compared to that of mCP. The sensitivity reach of FLArE with 10-tonne (red solid curve) and 100-tonne (orange solid curve) detector for each EM form factor is shown in Fig. 3. In the top panels, we present the results for the mass-dimension 5 operators, while the sensitivity reach for mass-dimension 6 couplings is shown in the bottom panels. In each case, current bounds are shown with gray-shaded regions. Light dark species χ with m χ O(10 MeV) are constrained dominantly by the past searches at the CHARM experiment [124,125] in the mass-dimension 5 case, and by supernovae bounds [69] for the mass-dimension 6 operators. Instead, heavier χs are constrained by mono-jet missing energy searches at the LHC [59,126]. In the plots, we also present other bounds on the considered effective operators from the BaBar [127], E613 [128], LSND [129], MiniBooNE [130], On the other hand, for mass-dimension 6 searches focusing on electron recoil are in general superseded by LEP and LHC missing energy searches. We also show the bound from a search for missing-energy in vector boson fusion (VBF) processes; however, it needs to be taken with caveat that, in this regime, the EFT based on Fµν -only is not the appropriate description of the underlying interaction [59].
We also show possibly more stringent constraints from di-jet searches at the LHC in which χs are produced via the vector boson fusion (VBF) process, cf. discussion in Ref. [59]. This production mode is enhanced in the presence of the χ couplings to only photons, although it also corresponds to the regime of this model in which unitarity might be violated. For this reason, we present this bound with a light gray color. We note that the relevant constraints become weaker than the aforementioned mono-jet bounds if the full hypercharge coupling is taken into account. In this case, however, further LEP bounds from invisible Z boson decays would again render the relevant region in the parameter space excluded.
In the case of mass-dimension 5 operators shown on top, we find that FLArE can probe comparable EM form factor couplings compared to the past proton-beam experiments such as CHARM-II, and yield competitive sensitivity with the LHC search based on the VBF production. In order to probe unconstrained regions in the parameter space of these models, a 100-tonne detector would be required. This will set new bounds for sub-GeV χ, while for larger χ masses the expected sensitivity becomes weaker as the dark states can no longer be effectively produced in rare J/ψ meson decays, cf. Fig. 1. In this case, missing energy searches at LEP and LHC become the most important. For reference, we also present the expected sensitivities of future proposed searches at the Belle-II, DUNE, and LDMX-II detectors, following Refs [57,65].
Instead, for the mass-dimension 6 operators, combination of SN1987A energy-loss bound [69] and constraints from missing energy searches at colliders [57,59] cover most of the parameter space down to a χ , b χ ∼ 10 −6 GeV −2 , while both FLArE and FLArE-100 only probe larger values of these couplings, a χ , b χ ∼ 10 −2 GeV −2 . In general, missing-energy searches do not require produced dark states to have recoil signals, thus they, naively, win over searches based on the scattering signals due to the lack of an additional squared coupling suppression, assuming the same intensity or luminosity.
We note, however, a relative difference between the mass-dimension 5 and 6 operators in the strength of these bounds compared to the FLArE sensitivity. This can be explained in a twofold way related to both the expected decreased FLArE χ detection prospects and increased sensitivity of missing-energy searches for χs coupled via the AM/CR couplings (dim-6) compared to lower-dimensional operators. As long as FLArE is concerned, this is first due to a suppressed χ-production rate which, in the AM/CR case, is dominated by heavier J/ψ mesons even for very forward-going χs, cf. right panels of Fig. 1 for comparison (note different normalization factors used for different couplings). Moreover, the E R -scaling of the differential scattering cross section for the mass-dimension 6 operators is the same as that for neutrino-electron interactions. As a result, in the search based on low-energy electron recoils, a good fraction of the BSM-induced signal is cut together with background events. Second, the prospects for missing-energy searches for χs coupled via the mass-dimension 6 operators are better than in the EDM/MDM case. This relies on the energy dependence of the χ production cross section for the process e − e + → χχγ which is increasingly proportional to the CM energy for higher-dimensional operators. These bounds will be further improved in the future HL-LHC searches [59].
We note that the FLArE detection prospects for χs coupled via the AM/CR operators could be improved by employing high-energy signal and nuclear scatterings, cf. Ref. [108] for a similar discussion for the light DM search in FLArE. Importantly, for these types of interactions, increased target mass enhances the χ scattering cross section, cf. Eq. (B2). While signal and background differentiation is less trivial in this case, in order to reduce systematic uncertainties, new physics events could be searched for as an excess in neutral current (NC) over charged current (CC) neutrino-induced events (NC/CC). We leave a detailed analysis of this effect for future studies. We do not expect, however, this search strategy to improve FLArE sensitivity by additional several orders of magnitude in the coupling strength which is needed to probe unconstrained regions in the parameter space of these models.
Last but not least, we note that χs coupled to photons via dimension 5 and 6 operators can also play the role of DM. In fact, they can be the thermal DM relic produced in the early Universe through the standard freeze-out mechanism; see solid black curves in Fig. 3 for the corresponding relic targets [57,59]. As can be seen, current bounds exclude the relevant regions of the parameter space of most of the considered scenarios for m χ < O(100 GeV) with a notable exception found for the MDM model. In this case, the FLArE and FLArE-100 detectors could probe values of the coupling constant for sub-GeV χs for which a subdominant DM thermal relic density is predicted. In the case of the EDM model, the currently allowed region in the parameter space to be probed in FLArE/FLArE-100 corresponds to too small values of the annihilation cross section and, therefore, predicts too large thermal DM relic density. In order to reconcile such scenarios with the cosmology, one would have to either modify the cosmological evolution of the Universe after χs freeze out or introduce changes in the model by, e.g., adding intermediate metastable particles into which χs could preferentially annihilate to suppress their abundance. This would likely introduce further bounds but possibly also detection prospects for such a model in the FPF which can go beyond the simplest scenario we focus on in this study.
We also stress that χ DM for mass-dimension 5 and 6 operators can also be probed in direct detection (DD) experiments and indirect detection (ID) searches, cf. Ref. [57,59] for recent discussion. Although we do not show this in Fig. 3, we note that these bounds can fully exclude the DM thermal freeze-out scenario of all the studied models for the mass range considered in this work. The only exception could be in the low-mass (sub-GeV) regime in which DD bounds become weaker, while ID constraints might be suppressed for a subdominant thermal DM relic density of χ. To go beyond the aforementioned constraints, one can consider the freezein mechanism that the χ relic is sourced from feeble interactions in the early universe [70].

VI. CONCLUSION AND DISCUSSION
The testing of the coupling strength between the dark sector species and SM photons plays a particularly important role in searches for new physics. This is first due to their possible impact on astrophysical and cosmological probes of such scenarios, but also because such couplings might naturally arise in UV models, albeit they are typically suppressed, consistent with current experimental bounds. It is then essential to study new ways to systematically constrain such couplings.
In this study, we have discussed the discovery potential of such searches in the far-forward region of the LHC. Specifically, we have focused on the recently proposed Forward Physics Facility [89][90][91] and the FLArE experiment [91,92] in which stable dark species could be searched for via their scatterings off electrons in the detector. We have studied the detection prospects in the search for popular milli-charged particles (mCPs), but also assuming that dark sector particles are electrically neutral and coupled to the electromagnetic field strength tensor F µν only via higher-dimensional operators. In the latter case, we have analyzed the EM form factors induced by mass-dimension 5 electric (EDM) and magnetic (MDM) dipole moments, as well as on the anapole moment (AM) and charge radius (CR) at mass dimension 6.
We find good detection prospects in the search for single-scattered electrons produced in interactions of the mCPs. In this case, FLArE can probe currently allowed regions of the parameter space of this model for the dark species mass in between 10 MeV and a few tens of GeV. We have also analyzed the mCP detection prospects in FLArE based on two coincident and collinear soft hits induced by mCP scatterings off electrons. This search strategy remains particularly important for the mCP mass of order tens of GeV or so. The FLArE searches could independently contribute to strengthen such expected bounds from the search based on low ionization signals in the FORMOSA experiment [12,91]. In addition, we stress that FLArE capabilities to study the statistics of signal events with different electron recoil energies could provide additional information to better reconstruct the mCP parameters in the case of their discovery during the HL-LHC phase. Together, this places the proposed FPF among the best-suited facilities to search for mCPs in a wide range of masses and, simultaneously, using various experimental techniques.
It is important to mention that FPF searches for mCPs in both FLArE and FORMOSA could gain from the presence of both detectors placed along the beam axis in the FPF. In particular, this could help with better triggering the events and with vetoing muon-induced backgrounds therefore further strengthening the expected future bounds. It also remains interesting to study potential FLArE discovery prospects for mCPs via observations of dim tracks in the LArTPC detector, cf. Ref. [40] for a similar discussion. We leave such analyses for future studies.
While the detection prospects in the FPF based on scattering signatures are very promising in models predicting low-energy depositions in the detectors, scenarios preferring large momentum transfer in the interaction are more difficult to disentangle from neutrinoinduced and other background sources. We illustrate this with the sensitivity reach of FLArE in the searches for EDM/MDM dark sector particles. This remains comparable to the current bounds from proton beam-dump experiments, while for the 100-tonne detector it can extend towards the future proposed searches in experimental facilities like Belle-II, DUNE, or LDMX-II. Even higher-dimensional operators leading to AM/CR interactions are more difficult to probe, with the expected FLArE reach not improving over the much stronger current bounds from the missing energy searches at LEP and LHC. In this case, an improved sensitivity reach could be obtained for more rich dark sector scenarios employing light new physics species and their couplings to the SM photons which could be studied via decay signatures in the FASER/FASER 2 experiments [79,81] and thanks to the secondary production processes, cf. Refs [48,134].
Last but not least, we stress that the results presented in the study have been obtained with the updated version of the FORESEE package [98]. On top of the previously implemented decay signatures, the package allows now for studying scattering signatures and contains further new physics models, as can be found in the public repository 2 . Besides the LHC, the package can also be used to analyze the detection prospects of mCPs and other dark sector species in potential far-forward searches at future hadron colliders with pp collisions at 27 where M is the mass of pseudoscalar meson, q 2 is the invariant mass square of χ-pair, θ is the angle between p χ in χ-pair rest frame and p χ + pχ in meson rest frame, and BR P →γγ is the branching ratio of pseudoscalar meson decaying into two photons. We adopt the values of BR P →γγ reported in PDG [139]. For vector meson decay, we implement the branching ratio into χ-pair for each millicharge and EM form factor rescaled from the branching ratio into electron-positron pair BR V →ee , which reads mCP: with M being the mass of vector meson and the benchmark values of BR V →ee taken from PDG [139].
Appendix B: DM-electron scattering cross section The differential scattering cross section for χe → χe reads with the functions g E (E R ) and g M (E R ) for all EM interactions can be found in [57]. Given that χ particles are highly-boosted for considered mass range, we can simplify Eq. (B1) assuming that E χ m e , m χ , E R and E R m e . For each EM form factor, the differential cross section that we adopt in FORESEE can be expressed as mCP: Assuming electrons are at rest in the lab frame, the maximal recoil energy E max R given E χ reads