Extending the Discovery Potential for Inelastic-Dipole Dark Matter with FASER

Neutral particles are notoriously difficult to observe through electromagnetic interactions. As a result, they naturally elude detection in most collider detectors. In this paper, we point out that neutral particles that interact through a dipole interaction can nevertheless be detected in far-forward detectors designed to search for long-lived particles (LLPs). In contrast to previous analyses that focused on neutral particles with elastic interactions, we consider inelastic interactions. This naturally leads to LLPs, and we demonstrate that FASER (and future experiments at the Forward Physics Facility) will be able to probe substantial regions of the associated parameter space. In particular, we find that FASER is capable of probing the region of parameter space wherein thermal freeze-out gives rise to an $\mathcal{O}$(GeV) dark-matter candidate with the appropriate relic abundance, as well as regions of parameter space that are difficult to probe at fixed-target experiments. FASER and its successor experiments may therefore play a critical role in the discovery of such a dark-matter candidate.

Since it began running almost 15 years ago, the Large Hadron Collider (LHC) has been one of the most effective tools we have for understanding the Standard Model (SM) and for probing physics beyond the Standard Model (BSM). Among its many successes, the LHC has discovered the Higgs boson and explored its properties in great detail. LHC data has also hinted at new physics in the flavor sector, and searches involving missing energy, displaced vertices, resonances, and a host of other phenomena have been used to constrain the properties of potential dark sectors [1][2][3]. However, despite this progress, the LHC has not yielded any definitive discovery of BSM physics.
Of course, new physics might be hiding within experimental blind spots. In particular, new physics might be found in the large-pseudorapidity regions, where there is a large flux of SM particles. The ForwArd Search ExpeRiment (FASER) [4][5][6] provides an ideal opportunity for probing this region. FASER is a cylindrical detector placed 476 m from the ATLAS interaction point (IP) along the beam axis. As such, FASER aims to capture the decay products of long-lived particles (LLPs). Shielded by 100 m of rock and concrete, FASER typically exploits the large SM activity near the IP as a means of dark-state production through meson decay. The FASER detector is already taking data, and an upgraded detector, known as FASER2, has been proposed, which would enhance the sensitivity of the current detector during the high-luminosity LHC (HL-LHC) era. This upgraded detector would be housed, along with several other experiments, within the proposed Forward Physics Facility (FPF) [7][8][9].
A common benchmark signature that FASER is well suited to probe is that of a dark photon decaying visibly to an e + e − pair. Such a signature can be used to constrain potential kinetic-mixing couplings of magnitudes ∼ 10 −5 − 10 −4 with dark-photon masses m A 1 GeV [4]. In addition to dark-photon models, FASER and the FPF have discovery prospects for a variety of new-physics scenarios in which an LLP, such as an axionlike particle, dark scalar, or heavy neutral lepton, decays to a pair of detectable SM particles. Although such scenarios are well studied in the literature (see, e.g., Ref. [5] and references therein), it is important to fully explore additional signatures to which FASER might also be sensitive.
In this paper, we examine the detection prospects for one such LLP signature at FASER -a monophoton signature that arises in scenarios in which the photon field couples to a pair of neutral BSM particles χ 0 and χ 1 with different masses m 1 > m 0 via an effective magnetic-dipole-moment (MDM) or electric-dipolemoment (EDM) operator. Such a scenario is particularly interesting from a phenomenological perspective, in part because it arises in the context of inelastic-darkmatter scenarios [10][11][12]. The inelastic MDM operator was also recently studied in the context of sterile neutrinos [13]. If χ 0 and χ 1 are sufficiently light, they can be produced together at a significant rate at the LHC via meson-decay processes involving such MDM and EDM operators. Moreover, these operators render χ 1 unstable, since χ 1 can decay via the process χ 1 → χ 0 γ. Thus, if the resulting χ 1 decay width is of the appropriate order, particles of this species emitted in the forward direction at the LHC could potentially decay inside the FASER detector. The calorimeter of the FASER detector can then detect the resulting monophoton signal [14,15].
From a theoretical perspective, inelastic-dark-matter models of this sort are also of interest because they represent the simplest realization of a non-minimal dark sector involving N dark states χ i with similar quantum numbers, where i = 0, 1, . . . , N − 1. The phenomenology of such dark sectors can differ significantly from that of minimal dark sectors, due to the presence of the additional dark-sector states (for a review, see Ref. [16]). Indeed, the presence of the additional dark-sector states generically gives rise to decay processes in which a heavier χ i particle decays into a final state involving one or more lighter χ i particles. Such decay processes alter the complementarity relations between different experimental probes of the dark sector [17,18]. Moreover, they can also lead to a variety of striking signatures at colliders [19,20] and even open up new ways of addressing the dark-matter problem, such as those that arise within the Dynamical Dark Matter (DDM) framework [21][22][23]. Non-minimal dark sectors involving multiple states with similar quantum numbers -and indeed often involving large values of N -emerge naturally in a variety of BSM contexts, including theories with extra spacetime dimensions, theories involving confining hidden-sector gauge groups, and string theory. A two-component inelastic-dark-matter model of the sort on which we shall focus in this paper may be viewed as the simplest example of a non-minimal dark sector in which these phenomenological possibilities arise.
One of the most important features of inelastic dark matter is that it provides a mechanism by which event rates at direct-detection experiments can be suppressed. Moreover, the dynamics involved in thermal freeze-out is significantly more complicated in such models, with annihilation, co-annihilation, up-and down-scattering, and decay processes all playing a role in generating the relic abundance of the lighter dark-sector state. Realizations of inelastic dark matter in which the dark-sector states couple to the visible sector via dipole-moment interactions in particular can give rise to a dark-matter abundance in accord with observation without running afoul of the applicable constraints [24][25][26]. Indeed, data from a variety of sources can constrain the parameter space of inelastic-dipole-dark-matter models, including data from collider experiments, fixed-target experiments, and gamma-ray telescopes. That said, there are also interesting regions of this parameter space in which existing searches have little or no sensitivity. Notably these include those in which the spectra are highly compressed and the lighter dark-sector state accounts for the entirety of the observed dark-matter abundance.
In this paper, we investigate the extent to which searches for a monophoton signature of χ 1 → χ 0 γ decays at FASER and FASER2 can constrain the parameter space of such inelastic-dipole-dark-matter models. As we shall see, these detectors provide coverage of this parameter space within certain regions that are not well constrained by existing searches. Indeed, FASER and FASER2 can probe parameter-space regions within which the splitting between the masses m 1 and m 0 of the two dark-sector states χ 1 and χ 0 is roughly ∆m ≡ m 1 −m 0 ∼ O(10 −2 )m 0 and where m 0 5 GeV. Remarkably, the reach of FASER and FASER2 within that parameter space includes regions wherein the relic abundance of χ 0 obtained from thermal freeze-out agrees with observed abundance of dark matter. We note also that the small mass-splittings ∆m that make the heavier state χ 1 longlived also suppress the energy of the photon into which it decays, thereby rendering the resulting monophoton signal challenging to detect at fixed-target experiments. At FASER, however, this suppression is compensated by the significant boosts provided to these decaying particles as a result of the high collision energies achieved at the LHC. Thus, as we shall see, FASER is capable of probing regions of parameter space that are difficult to probe at fixed-target experiments.
This paper is organized as follows. In Sect. II, we review the manner in which the dark-sector states couple to the visible sector in inelastic-dipole-dark-matter models. In Sect. III, we calculate the event rates for the signal and background processes relevant for monophoton searches at FASER and FASER2 within the context of these models and assess the discovery reach of these detectors within the parameter space of these models. In Sect. IV, we examine the annihilation and co-annihilation processes involved in the thermal freeze-out of χ 0 and χ 1 in the early universe and identify the regions of parameter space within which relic χ 0 particles constitute the majority of the dark matter at the present time. In Sect. V, we examine the additional considerations that constrain the parameter space of inelastic-dipole-darkmatter models in the region where FASER is sensitive. In Sect. VI, we conclude by summarizing our results and discussing possible directions for future work. Details of the calculation of meson branching fractions and thermal relic densities, along with a brief discussion of the detection prospects at SHiP [27,28], a proposed experiment complementary to FASER which will probe comparable lifetimes, are given in the Appendices.

II. INELASTIC-DIPOLE DARK MATTER
We consider an inelastic-dark-matter model in which the two dark-sector particles χ 0 and χ 1 are fermions that couple to the visible sector via an effective MDM or EDM operator in the interaction Lagrangian of the form where F µν is the field-strength tensor for the photon field and where Λ O , with O = {m, e}, denotes the corresponding operator-suppression scale 1 . We note that such effective operators can arise in a variety of UV-physics scenarios. For example, they can arise as a consequence of looplevel processes involving heavy particles with masses of O(TeV) or higher that carry SM charges. At scales above the electroweak-symmetry-breaking scale, these operators would presumably be replaced by operators with a similar structure, but with the U (1) Y and/or SU (2) L gauge bosons present in the unbroken phase of the SM in place of the photon [25]. However, this phase of the SM is not relevant for the phenomenological implications of inelastic-dipole dark matter that we consider in this paper. We note that these operators can also be viewed as arising from terms in the elastic tensor and pseudo-tensor currents J µν = χσ µν χ and J µν5 = χσ µν γ 5 χ involving a Dirac fermion χ that also acquires a Majorana mass [24][25][26]. Such a Majorana mass gives rise to a mass-splitting between the fermion mass-eigenstates -states that play the role of our χ 0 and χ 1 -and to operators with 1 An alternative parametrization is often used in the literature in which the strengths of the effective MDM and EDM interactions are parametrized by the dimensionful coupling coefficients µχ and dχ, which are related of our operator-suppression scales Λm and Λe by µχ = 2/Λm and dχ = 2/Λe. the inelastic coupling structures given in Eq. (2.1). We parametrize this mass-splitting in terms of the dimensionless ratio Since Dirac and Majorana terms in the dark-fermion mass matrix need not necessarily have a common origin, ∆ can in principle take a broad range of values. As such, we take the parameter space of our inelastic-dipole-darkmatter model to be characterized by the three parameters m 0 , ∆, and Λ O . The interactions in Eq. (2.1) have a variety of phenomenological implications. For example, in the presence of such interactions, χ 0 and χ 1 particles can be produced both at colliders and at beam-dump experiments via processes of the form ff → χ 0 χ 1 , where f denotes a SM fermion, via bremsstrahlung-like processes involving a virtual photon, or via the decays of vector and pseudoscalar mesons. Once produced, these particles can be detected in a variety of ways due to the up-scattering, down-scattering, and decay processes to which these same interactions give rise. Indeed, potential experimental signatures of inelastic-dipole dark matter include the scattering of χ 0 or χ 1 with atomic nuclei, events involving sizable missing transverse energy at colliders, and rare meson decays [29][30][31]. Existing limits on these processes place non-trivial constraints on our parameter space.
The decay process χ 1 → χ 0 γ provides an additional experimental probe of inelastic-dipole-dark-matter models. Indeed, searches for decays of this sort at B-factories and at the main LHC detectors likewise place non-trivial constraints on our parameter space [25]. However, as discussed in the Introduction, this same decay process can also potentially lead to signatures at FASER and FASER2 within regions of that parameter space that are largely unconstrained. One particularly interesting such region of parameter space turns out to be that in which ∆ ∼ O(0.01) and O(MeV) < m 0 < O(GeV). Within this region, not only can a significant flux of χ i pairs be generated in the forward direction, but the lifetime τ of χ 1 is also such that a significant fraction of χ 1 decays occur within the FASER detector. Indeed, as we shall justify below, the lab-frame decay lengthd = γβτ , where β and γ are the usual relativistic factors, can be expressed as a product of factors of the form where E χ1 is the energy of χ 1 . Given that FASER is located approximately 476 m away from the ATLAS detector, this expression ford makes it clear that the region of parameter space within which ∆ ∼ O(0.01) and O(MeV) < m 0 < O(GeV), while Λ O is near the TeV scale, is one in which the prospects of observing χ 1 decays insider FASER are particularly auspicious. Since the proposed location of FASER2 within the FPF lies a similar distance away from the ATLAS detector, this region of parameter space is auspicious for detection at FASER2 as well.
As we shall see, this region of inelastic-dipole-darkmatter parameter space is interesting not only because it is auspicious for detection at FASER; it is also interesting from a dark-matter perspective. Indeed, within this region the abundance of χ 0 obtained from thermal freeze-out and the subsequent decays of χ 1 can accord with the dark-matter abundance inferred from Planck data [32]. Moreover, since the value of the mass-splitting ∆m 10 keV within this parameter-space region of interest lies well above the lab-frame kinetic energy of a typical χ 0 particle of mass O(MeV) < m 0 < O(GeV) in the Milky-Way halo, event rates at direct-detection experiments are highly suppressed. As a result, such experiments are insensitive to dipole-interacting dark matter within this region of parameter space.

A. Production
A number of processes contribute to the overall production rate for χ 1 particles in the far-forward direction at the LHC in inelastic-dipole-dark-matter models. These include production from meson decay, Drell-Yan production, and secondary productioni.e., the production of χ 1 particles due to the up-scattering of χ 0 particles produced at the IP as they travel through the intervening material between the IP and the FASER detector [33].
The decays of neutral mesons produced at the AT-LAS IP turn out to provide the dominant contribution to this overall production rate. Indeed both vector and pseudoscalar mesons are produced at the LHC in copious amounts in the forward direction towards the FASER detector, which is coaxial with the ATLAS experiment. Moreover, the decays of both kinds of meson contribute to χ 1 production. Within the regime in which O(MeV) (2 + ∆)m 0 O(GeV), the relevant pseudoscalar mesons are π 0 , η, and η , while the relevant vector mesons are ρ, ω, φ, J/ψ, ψ(2S), and Υ(nS). The dominant contribution to χ 1 production from each of these meson species M is due to the decays of mesons that are produced in the far-forward directioni.e., with pseudorapidities η > η F , where η F = 9.2 and η F = 7.2 for FASER and FASER2, respectively -and have labframe energies E M ∼ O(TeV). Mesons with such energies are sufficiently highly boosted that the vast majority of χ 1 particles produced by their decays will likewise travel in the far-forward direction, regardless of the angle that the momentum vector of the χ 1 makes with the momentum vector of the decaying meson in the rest frame of the meson. Thus, we may obtain a rough estimate of the total cross-section σ (eff) χ1M for the production of χ 1 particles in the direction of FASER from the decay of a particular meson species M simply by multiplying the production cross-section σ M for mesons of that species with η ≥ η F by the overall branching fraction for M into final states including χ 1 . We obtain the σ M for each relevant meson species from the FORESEE code package [34].
A pseudoscalar meson P produces χ i particles primarily through the three-body decay process P → γγ * → χ 0 χ 1 γ. One could also consider the two-body decay P →χ 0 χ 1 , but the branching fraction for this process is loop-suppressed (because the leading diagram arises at one-loop order), helicity-suppressed, and further suppressed by Λ −4 O ∼ (10 TeV) −4 , so we shall not consider this production process further. By contrast, a vector meson V produces these particles primarily through the two-body decay process V → χ 0 χ 1 . The corresponding decay processes for the case of an elastic MDM or EDM interaction were considered in Refs. [29][30][31]. Full expressions for the branching fractions BR(V → χ 0 χ 1 ) and BR(P → χ 0 χ 1 γ) for the aforementioned vector-and pseudoscalar-meson decay processes are provided in Appendix A for both the MDM-and an EDM-interaction cases. We note that both BR(V → χ 0 χ 1 ) and BR(P → χ 0 χ 1 γ) are, to a good approximation, proportional to the square of the meson mass m M in the (2 + ∆)m 0 m M regime. Thus, although the cross-sections for the production of light pseudoscalars such as the π 0 in the forward direction at the LHC are significantly higher than those for much heavier vector mesons, the likelihood that each such vector meson will produce a χ 1 when it decays is enhanced simply by virtue of its being heavier. Moreover, the BR(P → χ 0 χ 1 γ) for all pseudoscalar mesons, regardless of their mass, experience an additional suppression relative to the BR(V → χ 0 χ 1 ) for vector mesons due to phase-space considerations. As a consequence of these two effects, we find that the decays of vector mesons in fact provide the dominant contribution to the production rate of χ 1 particles in the forward direction at the LHC within our parameter-space region of interest.
Given that vector-meson decays to χ i pairs dominate the production rate of χ 1 particles in the forward direction, it is instructive to compare the branching-fraction expressions BR MDM (V → χ 0 χ 1 ) and BR EDM (V → χ 0 χ 1 ) for these decays in the MDM-and EDMinteraction cases, respectively. In the regime in which m V (2 + ∆)m 0 , these two expressions coincide and are well approximated by the single expression where α is the fine-structure constant. By contrast, as the dark-fermion masses are increased to the point at which (2 + ∆)m 0 becomes comparable with m V , the expressions for BR MDM (V → χ 0 χ 1 ) and BR EDM (V → χ 0 χ 1 ) begin to deviate from each other appreciably. For example, in the elastic limiti.e., the limit in which ∆ → 0 -the ratio of these two branching-fraction expressions takes the particularly simple form (3.2) Thus, in this limit, the production rate of χ i pairs from vector-meson decay is larger for an MDM interaction than it is for an EDM interaction with the same value of Λ O . We find that this qualitative result likewise holds in the inelastic case in which ∆ = 0 as well. This can be understood as a consequence of angular-momentum conservation. In the decaying meson's center-of-mass frame, the MDM operator creates the χ 0 χ 1 final-state in the |S S z = |1 0 spin state, and since angular-momentum conservation requires J = 1, the final-state orbital angular momentum is L = 0. By contrast, the EDM operator creates a final state with spin state |S S z = |0 0 and L = 1. As a result, the partial width in the EDM case is P -wave suppressed, as can be seen in the squared amplitude in Eq. (A4), and thus FASER will have greater sensitivity for large m 0 for the MDM case.
In Fig. 1, we show the effective cross-sections σ (eff) χ1M for the production of χ 1 particles in the direction of FASER from the decay of several relevant meson species M as functions of m 0 for the choice of ∆ = 0.01. In order to display our results in as model-independent a manner as possible, we scale each effective cross-section by a factor of (Λ O /GeV) 2 to cancel the overall factor of Λ −2 O common to all of the σ (eff) χ1M . We observe that for this value of ∆, the decays of vector mesons collectively provide the dominant contribution to the production rate of χ 1 particles in the direction of FASER over the entire range of m 0 shown. Indeed, even for values of m 0 sufficiently small that the decay process π 0 → γχ 0 χ 1 is kinematically accessible, we find that the contribution to the total production rate from π 0 decay is subleading in comparison with the contributions from the decays of the heavier vector mesons ω and φ. Moreover, we note that these results are not particularly sensitive to the value of ∆, provided that ∆ 1. In addition to the contribution from meson decays, a contribution to the overall production rate of χ 1 particles at the LHC also arises from Drell-Yan processes of the form qq → χ 1 χ 0 . Since kinematic considerations imply that only mesons M with masses m M > (2 + ∆)m 0 contribute to the production of χ 1 particles through their decays, this Drell-Yan contribution can in principle be important in the regime in which (2 + ∆)m 0 is large. In order to determine the size of this latter contribution, we evaluate the differential cross-section d 2 σ/(dE χ1 dη χ1 ) for the Drell-Yan production of χ 1 particles as a function of their energy E χ1 and pseudorapidity η χ1 using the MG5 aMC@NLO [35] code package. We focus on the regime in which (2+∆)m 0 > 1 GeV, as parton-distribution functions are not well defined for momentum transfers below this scale. Based on the results of these simulations, we find that the contribution to the production rate of χ 1 particles from Drell-Yan processes is subleading in comparison with the contribution from meson decay.
Yet another contribution to the production rate for χ 1 arises due to secondary production [33]. Using MG5 aMC@NLO we evaluate this contribution by calculating the phase-space distribution of χ 1 after up-scattering via χ 0 + A → χ 1 + A off a nucleus A in the material located between the ATLAS IP and FASER. We incorporate appropriate form factors for the concrete, rock and for the FASERν tungsten target and calculate the event rate for subsequent χ 1 decay in FASER [33]. The region of parameter space for which this effect becomes potentially important is that within which 10 1 GeV Λ O 10 3 GeV and m 0 ∼ O(100) MeV. In particular we find that within this region of parameter space secondary production off tungsten in FASERν yields a potentially significant contribution to the signal rate. By contrast, we find that secondary production off the concrete and rock located between the ATLAS IP and FASER yields a subleading contribution to the signal rate.
While secondary production can indeed be important within this region of parameter space, this region turns out to be excluded by data from BaBar, CHARM-II, and LEP, as we shall discuss in detail in Sect. V. By contrast within regions of parameter space not excluded by current data, we find that the contribution to the total event rate is typically well below 1% of the contribution from meson decay.

B. Signal Rate
After a χ 1 particle is produced at ATLAS in our inelastic-dipole-dark-matter model, it decays via the process χ 1 → χ 0 γ. The lab-frame energy of the photon produced by the decay of a χ 1 particle with lab-frame energy E χ1 is where γ ≡ E χ1 /m 1 is the usual relativistic boost factor and where θ γχ1 is the angle in the rest frame of the χ 1 particle between the momentum vector of the photon and the axis along which the χ 1 particle is traveling. Since the χ 1 particle effectively inherits its boost factor from the decaying meson that produced it, and since the characteristic energy scale for the mesons is E M ∼ O(1 TeV), the characteristic size for this boost factor is γ ∼ E M /m M ∼ O(10 3 ). Thus, the characteristic photon energy for these decays within our parameter space is E γ ≈ γm 0 ∆ ∼ O(10 GeV). A single photon with such an energy can be detected by the calorimeter in either FASER or FASER2 [6].
The probability of observing a statistically significant number of monophoton events at FASER during Run 3 or at FASER2 during the HL-LHC era depends both on the geometry of the corresponding detector and on the decay properties of χ 1 itself. FASER and FASER2 each have (eff) χ 1 M for the production of χ1 particles in the direction of FASER from the decay of several relevant meson species M , shown as functions of m0 for ηF = 9.2. In order to display our results in as model-independent a manner as possible, we scale each cross-section by a factor of (ΛO/GeV) 2 to cancel the overall factor of Λ −2 O common to all of the σ (eff) χ 1 M . The solid curves correspond to the case of a MDM operator, while the dashed curves correspond to the case of an EDM operator. The results shown here correspond to the case in which ∆ = 0.01; however, we note that these results are not particularly sensitive to the value of ∆, provided that ∆ 1.
a cylindrical geometry in which the cylinder is coaxial with the beam-collision axis at the ATLAS IP. The components of FASER and FASER2 each include, starting from the side closest to the ATLAS detector, a dedicated decay volume, a tracking volume interleaved with a series of charged-particle trackers, a "pre-shower" station, and a calorimeter. Further details about each component are provided in Ref. [6]. The radius R of the cylindrical decay and tracking volumes, their combined length h along the axial direction of the cylinder, the distance L between the side of the decay volume closest to the ATLAS detector and the ATLAS IP, the value of η F that follows from the values of L and h, and the integrated luminosity L int expected at ATLAS during the corresponding running period are provided in Table I.
We note that the values of h given in Table I include both the dedicated decay volume of the detector and the tracking volume. Indeed, the detection of a single photon within the calorimeter of FASER or FASER2 does not rely on information from the tracker. Moreover, the scintillators separating the decay volume and tracking chambers of either detector are largely transparent to photons with E γ ∼ O(10 GeV) energies. Taken together, these two considerations imply that the tracking volume can also be viewed as an extension of the decay volume in searches for a monophoton signal of χ 1 → χ 0 γ decays within these detectors.
The characteristic decay length of χ 1 in the lab frame  is given byd where p χ1 and Γ χ1 denote the lab-frame momentum and proper decay width of χ 1 , respectively. Since all of the meson species that contribute significantly to the production rate of χ 1 particles decay promptly within the ATLAS detector,d represents the characteristic distance away from the IP that χ 1 particles travel before they decay. Within our parameter-space regime the proper decay width of χ 1 is well approximated by (3.5) in both the MDM and EDM cases. Since R h for both FASER and FASER2, we may approximate the distance that a χ 1 particle would traverse inside the decay volume of the detector if it didn't decay as h, regardless of its pseudorapidity η χ1 . Thus, the probability that a given χ 1 will decay inside FASER or FASER2 is approximately The total expected number of signal events from χ 1 decays within the FASER or FASER2 detector can therefore be approximated as where dσ χ1M /(dE χ1 dη χ1 ) is the differential cross-section for the production of χ 1 particles as a a function of E χ1 and η χ1 . We evaluate each of these differential crosssections numerically by integrating the product of the differential production cross-section for M obtained from the FORESEE package [34] and the appropriate phasespace factor over possible meson momenta.

C. Background
The primary detector backgrounds for an E γ ∼ O(GeV) monophoton signal of inelastic-dipole dark matter at FASER or FASER2 are those associated with muons or neutrinos produced at the ATLAS IP -or with muons produced by cosmic-ray showers -interacting with the calorimeter material [14,15,36].
Two scintillator stations located in front of the decay volume of FASER provide a veto on muons and other charged particles emanating from the direction of the AT-LAS IP, the efficiency of which is around 99.99% for each station individually [6]. As a result, the contribution to the overall background-event rate from such muons is expected to be negligible. An additional contribution to this overall rate also arises from events in which a muon is produced at the IP with a pseudorapidity η sufficiently large that it travels in the general direction of FASER, but sufficiently small that its path does not intersect either of the scintillator stations. If such a muon interacts with the external detector apparatus itself or some part of the surrounding infrastructure, this interaction can nevertheless produce a shower of particles. Depending on the directions in which these particles are emitted, some of these particles could potentially reach and deposit their energies in the calorimeter [37]. However, this contribution is also expected to be negligible, given that the particles produced in the shower are typically emitted in a narrow cone around the axis along which the muon was originally traveling.
The majority of cosmic-ray muons that pass through the FASER detector would likewise evade the scintillator veto. However, these muons also behave as minimum ionizing particles and typically give rise to energy deposits in the FASER calorimeter of O(MeV) -an energy scale far below the energy scale associated with the signal process. Moreover, the contribution to the backgroundevent rate from cosmic-ray muons can be further reduced by correlating timing information for calorimeter events at FASER with timing information for collision events at ATLAS. Given these considerations, this contribution is also expected to be negligible.
Finally, neutrinos produced in the far-forward direction at the ATLAS IP that undergo deep inelastic scattering with the material in the FASER calorimeter can give rise to energy deposits similar to those associated with the signal process. Background events of this sort can be distinguished from signal events on the basis of information from the high-granularity preshower detector [36,38]. Nevertheless, a residual contribution to the background-event rate still arises from two classes of neutrino events. The first class consists of events in which the scattering of the neutrino with the calorimeter generates a backsplash involving one or more photons that are subsequently detected by the pre-shower detector. We can significantly reduce the contribution to the background-event rate from the first class of events by imposing a maximum cut on the total energy deposited in the calorimeter, which is far higher -typically around O(100 GeV) [39] -for these neutrino-background events than it is for the events that constitute our signal. While we do not impose such a cut explicitly in our analysis, we note that the rate of these backsplash photon events is not expected to be high, even at energies above O(100GeV). Thus, we do not expect the imposition of such a cut to have a significant impact on our results.
The second class of background events consists of those in which a neutrino scatters in or near the final tracking layer and produces charged particles that are not detected by any of the trackers, but are nevertheless detected by the calorimeter. While a precise determination of this background would require a full detector simulation, the event rate for neutrino events of this class is expected to be extremely small. Nevertheless, a TeVneutrino scattering event is expected to produce O(10) charged particles with larger energies than our photon signal. Event-selection criteria designed to reject events which contain a significant number of such particles can therefore be effective in reducing this background.
Photons from χ 1 → χ 0 γ decay typically have energies E γ ∼ (1 TeV)∆ within our parameter-space region of interest. However, the FASER and FASER2 detectors are primarily designed to look for O(TeV) energy deposits, and the detection efficiency γ (E γ ) for a monophoton signal with E γ O(TeV) as a function of E γ is not well established. In order to illustrate the impact of this efficiency on the sensitivity of these detectors to our signal process, we focus for simplicity on the case in which γ (E γ ) can be modeled in terms of a binary detection threshold E γ,mini.e., we take within the (∆, E γ,min )-plane of the fraction f γ (E γ,min ) of signal events that would be detected at FASER for given E γ,min , relative to the number that would be obtained for E γ,min = 0. The solid, dashed, and dotted contours appearing in each panel of the figure correspond to different combinations of the parameters m 0 and Λ O . Two such contours -one corresponding to f γ (E γ,min ) = 0.1 and the other to f γ (E γ,min ) = 0.5 -are included for each parameter combination. The three panels of the figure simply display different regions of the (∆, E γ,min )-plane in order to illustrate the shape of these contours for the three m 0 and Λ O combinations shown therein. The results displayed in the figure correspond to the case of a MDM interaction. However, the value of f γ (E γ,min ) for an MDM interaction does not differ by more than 5% of the corresponding value obtained for an EDM interaction along the contours shown. In obtaining the results shown in Fig. 3, we have employed a cut on the photon energy of E γ > 300 MeV, which is well above expected energy deposits coming from minimum ionizing particles.

D. Results
In Fig. 3, we show contours of the integrated luminosity required to observe a monophoton signal of inelasticdipole dark matter at FASER (red curves) and FASER2 (blue curves) within the (Λ −1 O , m 0 )-plane for a few of different choices of ∆. Since the expected number of background events is effectively zero, as discussed in Sect. III C, we adopt the observation of at least three signal events as our discovery criterion. We emphasize that the luminosity contours can also be interpreted as event contours for a fixed total luminosity, L int (cf. Table I), i.e. for a contour, the number of events can be read as N χ1 = 3 × (L int /L) , where L is the luminosity of the contour. This can be useful in the case where there are unexpected backgrounds. The top, middle, and bottom panels in each column correspond to the choices ∆ = 0.001, ∆ = 0.01, and ∆ = 0.05, respectively. For smaller ∆, we expect to be sensitive only to values of Λ −1 O which lie within the region excluded by constraints from CHARM-II and LEP -constraints which are insensitive to ∆ for ∆ 1. While we do not consider smaller ∆ for this reason, we note that we expect FASER also to have some sensitivity for smaller ∆. The results in the left column correspond to the case of a MDM interaction, while the results in the right column correspond to the case of an EDM interaction. The dashed black curve in each panel indicates the contour along which the present-day relic abundance of χ 0 -a derivation of which shall be presented in Sect. IV -accords with the dark-matter abundance inferred from Planck data [32]. For ∆ 1, one expects our model to approach the elastic limit (i.e. ∆ = 0), and thus the relic target should remain fixed for ∆ < 10 −3 . The gray regions in each panel are excluded by existing constraints -constraints that will be discussed in more detail in Sect. V.
The overall shape of the sensitivity contours that appear in Fig. 3 can be understood as follows. For a given integrated luminosity and a given combination of m 0 and ∆, the range of Λ O for which a statistically significant number of χ 1 particles decay within the FASER or FASER2 detector is bounded from both above and below. The lower bound stems from the requirement thatd be sufficiently large that a non-negligible number of χ 1 particles reach the decay volume before they decay. Taken together, Eqs. (3.4) and (3.5) This implies that for a fixed value of ∆, we expect the contour in the (Λ −1 O , m 0 )-plane above which a significant number of χ 1 particles reach the FASER decay volume before they decay to take the form Λ −1 O ∝ m −2 0 . Indeed, the shape of the upper part of each sensitivity contour shown in each panel of Fig. 3 takes this form.
By contrast, the upper bound on the corresponding range of Λ O values for which FASER or FASER2 is sensitive for a given integrated luminosity and a given combination of m 0 and ∆ stems from the fact that for very large Λ O , most χ 1 particles that enter the decay volume of the detector pass through it without decaying. Ford L, the expression for the decay probability in Eq. (3.6) reduces to P (d, h, L) ≈ h/d. Moreover, the production rate of χ 1 , which is dominated by the contribution from vector-meson decay, is also suppressed for large As a result, the expected number of χ 1 decays within the FASER or FASER2 detector scales with Λ O , m 0 , and ∆ like Thus, for a fixed value of ∆, the maximum value of Λ O for which N χ1 drops below the sensitivity threshold scales with m 0 according to a proportionality relation of the form Λ O ∝ m 0 . This scaling behavior accounts for the shape of the lower part of each sensitivity contour shown in the figure.
The "bumps" that appear in these sensitivity contours -bumps which are the most prominent in the upper left panel -arise as a consequence of the decay kinematics of the individual meson species that contribute to the production of χ 1 . In particular, as m 0 increases, mesons with masses m M < (2 + ∆)m 0 become kinematically forbidden from decaying into χ i pairs, and as a result the overall production rate of χ 1 particles decreases.
Comparing the results shown in the corresponding panels in the left and right columns of the figure, we observe that the discovery reach obtained for an MDM operator structure at both FASER and FASER2 is comparable to, but nevertheless slightly larger than, the reach obtained for an EDM operator structure. This is a reflection of the fact that as m 0 approaches the kinematic cutoff m 0 = m V /(2 + ∆) above which the decay process V → χ 1 χ 0 is kinematically forbidden for a given vector-meson species, σ (eff) χ1M falls off more rapidly for an EDM operator structure than it does for a MDM operator structure, as indicated in Fig. 1. The primary message of Fig. 3, however, is that both FASER and FASER2 are capable of probing a sizable region of the parameter space of inelastic-dipole-darkmatter models that has not been meaningfully probed by other experiments. As one would expect, the discovery reach afforded by FASER2 within that parameter space is significantly larger than that afforded by FASER. However, we observe that both experiments are sensitive within hitherto unprobed regions of that parameter space wherein the χ 0 particle has the correct relic abundance to account for the entirety of the dark matter. Moreover, we emphasize that while regions of parameter space that lie below the black dashed line would overproduce dark matter in the context of the standard cosmology, such regions of parameter space would be viable in modified cosmological scenarios involving, for example, additional entropy production from the out-of-equilibrium decays of heavy particles after thermal freeze-out, but before the nucleosynthesis epoch.

IV. RELIC ABUNDANCE
In order for the χ i to play the role of the dark matter in our inelastic-dipole scenario, the collective present-day abundance of these states must be consistent with the dark-matter abundance Ω DM ≈ 0.26 inferred from Planck data [32]. Thermal freeze-out provides a natural mechanism for generating abundances for both χ 0 and χ 1 in this scenario. Within the region of parameter space relevant for FASER, χ 1 is sufficiently short-lived that essentially all χ 1 particles present in the universe at the end of the freeze-out epoch decay to χ 0 particles well before present time -and for that matter, well before the Big-Bang Nucleosynthesis (BBN) epoch. Thus, within this region of parameter space, only χ 0 contributes to the presentday abundance of dark matter. Nevertheless, as we shall see, co-annihilation processes involving χ 1 particles during the freeze-out epoch can have a significant impact on the resulting χ 0 abundance at late times.
The evolution of the number densities n i of the χ i in the early universe is governed by set of coupled Boltzmann equations, each of which takes the form where H is the Hubble parameter and where C i represents the contribution to the rate of change of n i from scattering and decay processes. Equivalently, these differential equations may be formulated in terms of the comoving number densities Y i ≡ n i /s of the dark-sector species, where s denotes the entropy density of the universe. Indeed, during periods wherein the total entropy S ∝ sa 3 within a comoving volume is conserved, the Boltzmann equations for the Y i each take the particularly simple form In inelastic-dipole-dark-matter models, the processes that can contribute significantly to the C i include the following, where = e, µ, τ denotes a charged-lepton species: • t-channel annihilation processes χ i χ i ↔ γγ; • s-channel co-annihilation processes χ i χ j ↔ + − and χ i χ j ↔ hadrons, where i = j; • up-and down-scattering processes χ 0 ± ↔ χ 1 ± and analogous processes involving hadrons; • decay and inverse-decay processes χ 1 ↔ χ 0 γ.
The annihilation and co-annihilation processes listed above can change the total comoving number density Y tot ≡ Y 0 + Y 1 of dark-sector particles. By contrast, the other processes listed above can transfer comoving number density from one dark-sector species to the other, but have no effect on Y tot .
We can express the C i for each χ i as a sum of four terms, each associated with one of the four classes of process itemized above. In particular, C 0 and C 1 may be written in the form where n eq i is the equilibrium number density of χ i , where n eq is the equilibrium number density of the chargedlepton species , and where Γ denotes the thermallyaveraged decay rate for χ 1 . In addition, σ 00 v , σ 11 v , σ 12 v , and σ 0 v denote the thermally-averaged crosssections for the annihilation of two χ 0 particles, the annihilation of two χ 1 particles, the co-annihilation of a χ 0 and a χ 1 particle, and the up-scattering process χ 0 ± → χ 1 ± , respectively. Expressions for the thermally-averaged cross-sections and decay widths for all relevant processes are provided in Appendix B.
Within our parameter-space region of interest, χ 0 and χ 1 always depart from thermal equilibrium only after they have both become non-relativistic. Within this region of interest, whether or not co-annihilation processes have a significant impact on the overall dark-matter abundance depends primarily on the value of ∆. When ∆ 1, co-annihilation processes are relatively unimportant because the number density of χ 1 will be exponentially suppressed relative to that of χ 0 at the time when χ 1 effectively becomes non-relativistic. As a result, the cosmological evolution of χ 1 is essentially decoupled from that of χ 0 . By contrast, when ∆ 1, the number densities of χ 0 and χ 1 remain comparable until much later, and the effect of co-annihilation processes cannot be ignored. Indeed, within the region of parameter space relevant for FASER, wherein ∆ 1, co-annihilation processes turn out to play a significant role -and indeed often a dominant role -in determining the overall present-day abundance of χ 0 .
Of course, the freeze-out dynamics in this scenario, and therefore the late-time abundance of χ 0 , depends not only on ∆m, but on m 0 as well. Since the temperature at which thermal freeze-out occurs is typically m 0 /O(10), for small masses m i m µ the production of χ i involves only electrons and positrons. However, for m i m µ , processes involving muons are also relevant. At even higher temperatures, interactions involving hadrons or tau leptons become relevant as well.
At first glance, it might seem that all of the collision terms in Eq. (4.3) contribute to maintaining thermal equilibrium between the χ i and the thermal bath of SM particles, given that all of these terms vanish when n i = n eq i . Were this the case, the up-and downscattering process in which the χ i participate, as well as the decay and inverse-decay processes, would serve to keep Y 0 and Y 1 in equilibrium for a significant duration, given that the number densities of the light SM particles that participate in these processes are still quite large during the freeze-out epoch. As a result, the relic abundance of χ 0 would be heavily suppressed. However, as we have discussed above, these processes -unlike annihilation and co-annihilation processes -have no effect on Y tot . Thus, after the annihilation and co-annihilation rates drop below H, both χ 0 and χ 1 depart from the thermal equilibrium with the SM thermal bath. Nevertheless, the up-and down-scattering processes, together with the decay and inverse-decay processes, serve to maintain the relations between Y 0 and Y 1 even after this departure from equilibrium occurs. These relations imply that Thus, at late time, when T m 1 − m 0 and the equilibrium abundances therefore satisfy Y eq 1 Y eq 0 , it therefore follows that Physically, the results in Eq. (4.6) reflect the fact that even if Y 0 and Y 1 are comparable at the moment when χ 0 and χ 1 initially depart from thermal equilibrium with the radiation bath, up-scattering and inverse-decay processes will become increasingly kinematically disfavored as T drops. As a result, the initial population of χ 1 particles present at the point at which this departure occurs will gradually be converted into χ 0 particles via decay and down-scattering processes until the majority of this initial population is all but depleted and χ 0 freezes out.
In Fig. 4, we show how the comoving number densities Y 0 and Y 1 evolve as a function of m 0 /T for an example choice of model parameters that yields a relic abundance for χ 0 that accords with Ω DM . Since ∆ 1 for this choice of parameters, the equilibrium comoving number densities Y eq 0 and Y eq 1 (indicated by the solid orange and dashed cyan curves in the figure) are very similar. As a result, the actual comoving number densities Y 0 and  and Y eq 1 . The dotted gray line indicates the value of Yi that yields the correct present-day dark-matter abundance for this choice of model parameters. The yellow curve represents the quantity Y0(Y eq 1 /Y eq 0 ), which, in accord with Eq. (4.4), effectively coincides with Y1 across the entire range of m0/T shown. Y 1 , which track their equilibrium values very closely at high temperatures, likewise remain quite similar at high temperatures. However, at around m 0 /T 20, when Y 0 and Y 1 begin to depart appreciably from their equilibrium values, they also begin to differ appreciably from each other. Indeed, by the time at which Y 0 has effectively settled into its late-time asymptotic value -a time that corresponds roughly to the point at which m 0 /T ∼ 200 -this comoving number density is already an order of magnitude larger than Y 1 . We also observe that Y 1 tracks the quantity Y 0 (Y eq 1 /Y eq 0 ), in accord with the result in Eq. (4.4).
We also note that the value of Λ O required to obtain the correct thermal relic density in the EDM case is smaller than it is in the MDM case for the same choice of m 0 and ∆. This is ultimately a consequence of angularmomentum conservation. In the regime in which ∆ 1 and χ 0 and χ 1 are nearly degenerate in mass, the thermal relic density is determined primarily by co-annihilation processes such as χ 0 χ 1 → e + e − , which are mediated by an s-channel photon. In the MDM case, these are Swave processes: as noted in Sect. III A, the initial state has spin S = 1, and since the final state also has spin S = 1, angular-momentum conservation can be satisfied by states with orbital angular momentum L = 0. By contrast, in the EDM case, they are P -wave processes: the initial state has spin S = 0 and the final state has spin S = 1, so angular-momentum conservation requires orbital angular momentum L = 1. The value of Λ O must therefore be smaller in the EDM case to compensate for the resulting P -wave suppression.
We can obtain a more quantitative picture of how this P -wave suppression affects the freeze-out dynamics in the EDM case relative to the MDM case by comparing the corresponding expressions for the cross-section σ MDM 01→ee and σ EDM 01→ee for the co-annihilation process χ 0 χ 1 → e + e − -expressions respectively obtained by substituting the C (s) factors in Eqs. (B9) and (B9) into Eq. (B8). In particular, since ∆ 1 within our parameter-space region of interest, one can gain some insight into how the differences between these expression affect our results by considering how these expressions behave in the ∆ → 0 limit. In this limit, one finds that these expressions reduce to where v = 1 − 4m 2 0 /s is the speed of either initialstate particle in the center-of-mass frame. Thus, we see that the co-annihilation cross-section in the EDM case is suppressed relative to the co-annihilation cross-section in the MDM case by a factor of v 2 . At freeze-out, when v 2 ∼ 0.1, the two cross-sections are equal when Λ m ∼ 5Λ e . Indeed, this is reflected in the locations of the Ω χ0 = Ω DM contours appearing in Fig. 3.

V. EXISTING CONSTRAINTS
The parameter space of inelastic-dipole-dark-matter models is constrained by a variety of experimental and observational considerations. Within the region of that parameter space that we have identified as being the most relevant for FASER, the most important such constraints are those from colliders and from fixed-target experiments. In addition, in cases in which the present-day abundance of χ 0 particles is non-negligible, astrophysical considerations also imply constraints on that region of parameter space. We discuss these constraints in turn below.

A. Decay Signatures: Colliders
The BaBar experiment at the PEP-II e + e − facility is sensitive to regions of parameter space near the thermalrelic contour for inelastic-dipole dark matter and indeed can be seen as complementary to FASER. The resulting constraints on this parameter space were investigated in Ref. [25] for the case of a MDM interaction. We extend this analysis here both to the case of an EDM interaction and to smaller values of ∆.
Limits on events involving either one or two energetic photons and sizable missing energy / E at BaBar place constraints on inelastic-dipole dark matter. The darksector particles χ 0 and χ 1 can be produced there either through the t-channel process e + e − → γχ 1 χ 0 or through the s-channel process e + e − → χ 1 χ 0 . The leading contribution to the signal-event rate in the 2γ + / E channel arises from t-channel production, followed by the prompt decay of the χ 1 . By contrast, sizable contributions to the signal-event rate in the γ + / E channel arise both from schannel production, followed by the prompt decay of the χ 1 , and from t-channel production, followed by the χ 1 escaping the detector before it decays.
For our analysis of the signal rate in both the 2γ + / E and γ + / E channels, we make use of the BaBar monophoton trigger, which was in effect while 60 fb −1 of integrated luminosity was accumulated. This trigger requires that the leading photon in the event have energy E γ ≥ 2 GeV. For both channels, we also require that / E ≥ 20 MeV. For the 2γ + / E channel, following Refs. [25,29], we also require that an additional photon with E γ ≥ 50 MeV be present in the event.
For the 2γ + / E channel, we assume no background and require that the χ 1 decays promptly (i.e., within 1 cm of the interaction point). Since any χ 1 particles which might have been detected at BaBar would necessarily have been produced with much smaller boosts, this prompt decay length is the one relevant for our parameter-space region of interest. By contrast, sizable SM backgrounds to the γ + / E signature arise from events involving neutrino production via neutral-current interactions and from events in which a photon is produced in conjunction with additional photons or e + e − pairs that escape detection. We find that the constraints derived from the γ + / E channel turn out to be subdominant to the constraints from the 2γ + / E channel, and thus we focus on these constraints in what follows.
The signal contribution in the 2γ + / E channel involves the prompt decay of the χ 1 , and so the corresponding constraints on our model-parameter space extend to arbitrarily small Λ O . These constraints exclude the shaded region of model-parameter space labeled "BaBar" in each panel of Fig. 3. We note that the Belle-II experiment [40] at SuperKEKB, which is currently taking data, will be able to extend the reach of B-factory experiments within this parameter space to larger Λ O in the near future.
Since χ 1 can also be produced via analogous processes at hadron colliders, LHC searches likewise place constraints on the parameter space of inelastic-dipole-darkmatter models. The extent to which such searches are capable of probing the parameter space of such models was investigated in Ref. [25] for the case of a MDM interaction. While the results of this analysis indicate that LHC searches are in general capable of probing that parameter space down to Λ O ∼ O(10 3 ) GeV for a broad range of m 0 and ∆, they also indicate that for ∆ 0.05, the region of parameter space that these searches are capable of probing does not extend beyond the region that is already excluded by BaBar data. Thus, the constraints from LHC searches on the region of that parameter space relevant for FASER are subdominant to those from BaBar and need not be considered further.
One could also consider constraints from the LEP experiment constraining invisible Z decays. Such constraints were considered in Ref. [25], where the authors explored a similar inelastic dipole operator, but with a coupling to the hypercharge gauge boson instead of the photon. For the magnetic dipole case, they derived a lower bound Λ m 10 3 GeV on the effective scale associated with the corresponding operator. Depending on the UV theory, a similar constraint might be applicable to our model as well. However since the χ i have no direct coupling to the physical Z boson in the effective theory we consider in this paper, we do not consider this constraint further.

B. Decay Signatures: Fixed-Target Experiments
A variety of experiments in which a beam of SM particles -typically electrons or protons -collides with a fixed target located some distance in front of a particle detector impose non-trivial constraints on inelasticdipole-dark-matter models. Indeed, χ 1 particles can be produced by these collisions, and some fraction of these χ 1 particles subsequently decay within the detector, producing photons that can then be detected. At electron beam-dump experiments, χ 1 particles are produced predominately via bremsstrahlung-like processes such as e − N → e − N χ 1χ0 , where N denotes an atomic nucleus. At proton beam dumps, these particles are produced primarily via the decay of various mesons produced by the collisions between the incident protons and the nuclei in the target. We shall consider the constraints from these two classes of fixed-target experiments in turn.
Among electron beam-dump experiments, E137 [41] and E141 [42] are the most relevant for constraining the region of model-parameter space pertinent for FASER. At each of these experiments, electrons with lab-frame energies E e ∼ O(10) GeV collided with a dense target -aluminum in the case of E137, tungsten in the case of E141. These experiments also had similar, O(GeV) detection thresholds. By contrast, the beam-dump experiment at Orsay [43], which had a similar detection threshold but a far lower electron-beam energy E e = 1.6 GeV, yields far weaker constraints within this same region of parameter space. The beam-dump experiment at KEK [44] had an electron-beam energy E e = 2.5 GeV, which is also much lower than that of either E137 or E141. However, this experiment also had a much lower detection threshold -around O(100 MeV). Nevertheless, given the significant background from photon bremsstrahlung at energies E γ 1 GeV, this lower detection threshold is unlikely to result in a significant improvement in sensitivity to monophoton signals. Moreover, the total number of electron collisions attained over the lifetime of the KEK experiment is a factor of 10 3 smaller than the number attained at E137. As a re-sult, we do not expect the KEK experiment to yield constraints competitive with those from E137 within our parameter-space region of interest.
We assess the sensitivity of E137 and E141 within the parameter space of our inelastic-dark-matter model in a manner analogous to that employed in similar studies of other LLP scenarios in Refs. [42,45,46]. In particular, we simulate the production of photons at these experiments using the MG5 aMC@NLO code package [35] and incorporate the appropriate form factors for the target material in each experiment -form factors that account for the effects of coherent scattering as well as nuclear and atomic screening [46][47][48]. We find that neither E137 nor E141 is capable of achieving a statistically significant number of events necessary to constrain the region of model-parameter space relevant for FASER, even for zero background. This is primarily due to the fact that kinematic considerations place an upper bound E γ ≤ [E e + (E 2 e − m 2 0 ) 1/2 ]∆ on the lab-frame energy of a photon produced by χ 1 → χ 0 γ decay at fixed-target experiments of this sort. Since ∆ 0.1 within our parameter-space region of interest, E γ typically lies significantly below the threshold for detection. Indeed, it is only for ∆ > 0.1 that the E137 and E141 data become constraining for inelastic-dipole dark matter.
Among proton beam-dump experiments, Nu-Cal [49,50] and CHARM-II [51,52] are the most relevant within our parameter-space region of interest. The Nu-Cal experiment, which primarily functioned as a neutrino detector, collided approximately 10 18 protons -each with a lab-frame energy of E p ≈ 69 GeV -with an iron target over the course of its operation. Scalar and vector mesons produced in these collisions could decay into χ 1 particles via the processes discussed in Sect. III A. Decays of these χ 1 particles to photons within the Nu-Cal detector could then give rise to a signal.
In assessing the expected number of signal events at Nu-Cal in inelastic-dipole-dark-matter models, we adopt the procedure for modeling the scalar-and vector-meson production cross-sections used in Ref. [50]. We relate the inclusive cross-section σ pFe → π 0 X for π 0 production via proton scattering with an iron nucleus N to the inclusive cross-section for π 0 production via proton-proton scattering where X denotes any additional SM particles in the final state, where A =56 is the atomic mass number of iron, and where α = 0.55. We also model the differential crosssection for π 0 production via proton-proton scattering as an average of the corresponding differential cross-sections for π + and π − production, which have been measured empirically [53]. In other words, we take where z is the fraction of the outgoing momentum of the incident proton carried by the π 0 and where p T is the magnitude of the component of its momentum vector transverse to the beam axis. Finally, following Ref. [31], we estimate the differential cross-sections for the production all heavier neutral scalar-and vector-meson species M as where the scaling factor N values for all relevant meson species were computed in Ref. [31] using Pythia 8.2 [54,55] simulations, and we adopt these N (POT) M values as well. Given these differential cross-sections, we may calculate the spectrum of χ 1 particles produced in the direction of the Nu-Cal detector by the decays of these mesons in a straightforward manner. This detector consisted of a cylindrical decay volume 23 m in length and 2.3 m in diameter located 64 m behind the target with its axis aligned with the axis of the proton beam. In assessing the total expected number of signal events from χ 1 → χ 0 γ decays at Nu-Cal for any combination of Λ O , m 0 , and ∆, we require not only that the χ 1 decay within this decay volume, but also that the lab-frame energy of the resulting photon satisfy E γ > 3 GeV -a criterion that derives from the fact that Nu-Cal was only capable of distinguishing an electromagnetic from a hadronic shower when the total energy associated with the shower exceeds this threshold.
By comparing the expected number of signal events for different combinations of Λ O , m 0 , and ∆ to the number of events satisfying the same criteria which were actually observed at Nu-Cal during its full run -a total of 5 such events were observed, and the expected background was 3.5 events -we can derive constraints on our modelparameter space. These constraints exclude the shaded region of model-parameter space labeled "ν-Cal" in the bottom two panels of Fig. 3, which correspond to the parameter choice ∆ = 0.05. Similar exclusion regions do not appear in the other panels of the figure, however, because too few photons satisfy the E γ > 3 GeV criterion for such small values of ∆. Since the vector mesons ρ and ω, which have the highest production rates at Nu-Cal, are kinematically forbidden from decaying to χ i pairs for m 0 390 MeV, the results of this experiment have little constraining power within our parameter space for m 0 above this value.
The CHARM experiment and its successor CHARM-II [51,52] would likewise have been able to produce χ i pairs through meson decay. Thus, the results of these experiments also constrain the parameter space of inelasticdipole-dark-matter models. The original CHARM experiment collided approximately 7.1 × 10 18 protons -each with an energy of E p = 450 GeV -with a copper target over the course of its operation. The center of the detector, which was situated 487.3 m behind the target, had a transverse area of 2.4 × 2.4 m 2 . Following [56], which searched for a monophoton signal, we take the fiducial decay volume to be 2.4 × 2.4 × 11.2 m 3 . The CHARM II experiment collided approximately 2.5 × 10 19 protonseach also with an energy of E p = 450 GeV -with a beryllium target over the course of its operation. The detector, which was situated 870 m behind the target, had a transverse area of 3.7 × 3.7 m 2 . We take the length of the fiducial decay volume to be 35.7 m, which represents the length of the entire calorimeter, [57] in the direction parallel to the beam axis. While the distance between the target and the detector at each of these experiments was similar to the distance between FASER and the ATLAS IP, the characteristic lab-frame energy of a χ 1 particle produced by collisions of the proton beams at CHARM and CHARM-II with their respective fixed targets was significantly lower than the characteristic lab-frame energy of a χ 1 particle produced in the forward direction at the LHC. As a result, CHARM and CHARM-II data constrain larger values of Λ O than does FASER data, which lead to shorter proper lifetimes for χ 1 .
In deriving the constraints from CHARM data on our parameter space, we largely follow the procedure outlined in Refs. [31,56]. Likewise, in assessing the corresponding constraints from CHARM-II data, we largely follow the procedure outlined in Ref. [31]. For each experiment, we derive a differential cross-section for π 0 production of the Bonesini-Marchionni-Pietropaolo-Tabarellide-Fatis (BMPT) form [58] using the BdNMC code package [59]. We then obtain the differential cross-sections for the production of η, η , ρ, ω, φ and J/ψ by scaling this differential cross-section for π 0 production by an overall factor N (POT) M , as in Eq. (5.3). From these differential cross-sections, we may derive an expected number of signal events from χ 1 decays within the decay volume at either CHARM or CHARM-II for a given combination of Λ O , m 0 , and ∆ in much the same way as we derived the corresponding signal-event count at Nu-Cal. In accord with the cuts performed on the energy of the electromagnetic shower in Ref. [60], we require that the lab-frame energy of the photon satisfy 7.5 GeV ≤ E γ ≤ 50 GeV for the CHARM analysis. Likewise, for the CHARM-II analysis, we require that 3 GeV ≤ E γ ≤ 24 GeV in accord with the corresponding cuts performed in Ref. [51,57]. Although the minimum value of E γ for both of these windows is considerably higher than the threshold value employed in the Nu-Cal analysis above, E p is also considerably higher for both CHARM and CHARM-II than it is for Nu-Cal. As a result, the characteristic energies of the χ 1 particles are also significantly higher, and a non-negligible fraction of photons produced by χ 1 → χ 0 γ decay survive these E γ cuts, even for ∆ = 0.01.
We derive constraints on our model-parameter space by comparing the expected number of signal events obtained for different combinations of Λ O , m 0 , and ∆ to the number of events that were actually observed at CHARM or CHARM-II. However, there are significant uncertainties in the expected backgrounds at these experiments [51,56,60]. Thus, since our primary aim in this paper is to identify regions of our parameter space within which FASER and FASER2 are sensitive, but which are clearly not constrained by other experiments, we maximize the size of the exclusion regions from CHARM and CHARM-II by assuming a background of zero events at each experiment. The resulting constraint from CHARM-II data excludes the shaded region of model-parameter space labeled "CHARM-II" in each panel of Fig. 3. Since we find that these regions completely subsume those excluded by CHARM data, we do not include separate exclusion contours for CHARM in this figure.

C. Constraints from the Elastic Limit
In addition to the signals discussed above, which arise as a consequence of χ 1 decay and therefore vanish in the elastic limit, i.e., the limit in which ∆ → 0, there are a number of potential signatures of dipole-interacting dark matter at colliders, fixed-target experiments, and neutrino facilities that are associated with processes that do not vanish in that limit. These signatures include, for example, electron or nucleon recoils precipitated by collisions of the χ i with the detector material in fixed-target experiments and monojet + / E signatures at hadron colliders. The non-observation of such signatures at experiments including E613 [61], MiniBooNE [62], LSND [63], LEP [64], CHARM-II [51], and LSND [63] places nontrivial constraints on inelastic-dipole-dark-matter models.
Since the mass-splittings m 1 − m 0 within our parameter-space region of interest are quite small in comparison with the characteristic energy scales of these experiments, the corresponding constraints are well approximated by those obtained for a purely elastic dipole interaction. Detailed studies of the experimental constraints on elastic dipole-interacting dark matter have been performed for both the MDM and EDM cases [29][30][31]. Throughout most of our parameter-space region of interest, the leading such constraints are those from limits on electron recoils at CHARM-II or from limits on monophoton + / E searches at LEP. The shaded region of model-parameter space labeled "CHARM-II + LEP (∆-Ins)" in each panel of Fig. 3 is excluded by constraints of this sort. These constraints collectively impose a bound Λ O 1 TeV within our parameter-space region of in-terest that is not only independent of ∆, but also not particularly sensitive to the value of m 0 .

D. Astrophysical Constraints
In cases in which χ 0 particles have a significant presentday abundance, their annihilation within the halos of galaxies via the t-channel process χ 0 χ 0 → γγ can give rise to a monochromatic photon signal at gamma-ray telescopes. The leading constraints on such a signal are currently those from Fermi-LAT data, which imply an upper bound on the corresponding thermallyaveraged annihilation cross-section that varies with m 0 from σv 10 −34 cm 3 s −1 for m 0 ≈ 1 MeV to σv 10 −29 cm 3 s −1 for m 0 ≈ 1 GeV [65]. This thermallyaveraged cross-section may be estimated as (5.4) which is well below the Fermi-LAT bound for all m 0 within our parameter-space region of interest. Constraints from other instruments (such as EGRET) on monochromatic gamma-ray signals from the annihilation of sub-GeV dark-matter particles have also been assessed [66], but turn out to be subleading in comparison with the Fermi-LAT bounds.
In principle, an additional constraint on inelasticdipole-dark-matter models arises from observational limits on distortions of the cosmic microwave background [67] due to annihilation, co-annihilation, or scattering processes involving the χ i during or after recombination. However, we find that these limits do not constrain any portion of our parameter-space region of interest that is not excluded by other constraints. Indeed, since the cross-section for the annihilation process χ 0χ0 → γγ is proportional to Λ −4 O , the corresponding annihilation rate is highly suppressed for values Λ O 1 TeV that are consistent with constraints from collider and fixed-target experiments. Moreover, the lifetime of the χ 1 is sufficiently short throughout our parameter-space region of interest that the abundance of χ 1 particles is negligible at the time of recombination. At the same time, for m 0 1 MeV and ∆ 0.001, the mass-splitting between χ 0 and χ 1 satisfies m 1 − m 0 1 keV. Since this significantly exceeds the temperature T ∼ O(10 eV) during the recombination epoch, the production of χ 1 by up-scattering processes immediately prior to or during recombination is highly Boltzmann-suppressed. Taken together, these considerations imply that the co-annihilation rate during recombination is negligible.
Limits on the effective number N eff of neutrino species during the nucleosynthesis epoch place stringent constraints [68] on light thermal relics with highly suppressed couplings to the visible sector. However, such constraints only become relevant for m 0 10 MeV, and do not im-pact the region of model-parameter space relevant for FASER.
Finally, direct-detection experiments also place constraints on the scattering cross-sections of χ 0 with atomic nuclei or with electrons. However, inelastic-scattering rates for a dark-matter particles with m 0 ∼ O(GeV) at such experiments are highly suppressed for ∆ 10 −6 . While elastic scattering via loop-level processes still yield a contribution to the overall scattering rate for larger values of ∆ [69], this contribution is far too small within our parameter-space region of interest to be constraining.

VI. CONCLUSIONS AND DIRECTIONS FOR FUTURE RESEARCH
In this paper, we investigated the prospects for detecting a signature of inelastic-dipole dark matter at the FASER detector and its proposed upgrade FASER2. This signature involves the production of the heavier of the two dark-sector states present in the theory by collisions at the ATLAS IP and its subsequent decay into the lighter dark-sector state and a single photon inside the FASER decay volume. We determined the discovery reach of FASER and FASER2 within the parameter space of this model for the case of either a magnetic or an electric dipole-moment interaction and showed that these detectors are capable of probing regions of this parameter space within which other experimental probes have little or no sensitivity. Moreover, these include regions wherein thermal freeze-out yields an abundance for the lighter dark-sector state that agrees with the observed abundance of dark matter.
Several comments are in order. First, it is worth noting that one of the reasons why FASER and FASER2 are capable of probing regions of inelastic-dipole-dark-matter parameter space within which existing experiments have little constraining power is simply the far higher energy of the LHC. Indeed, as we saw in Sect. V, one of the reasons why current fixed-target bounds are not terribly constraining within our parameter-space region of interest is that the majority of photons produced by χ 1 decays within the detector at these experiments would have had lab-frame energies that are below the corresponding detection thresholds. By contrast, the substantial boost provided to χ 1 particles at the LHC ensures that the photons produced by the decays of these particles within the FASER detector have far higher energies. Moreover, the time-dilation factor provided by these large boosts enhances the ability of FASER to probe regions of parameter space in which the lifetime of χ 1 is considerably shorter. Since small mass splittings are a generic possibility that leads to long-lived particles, and since these small mass splittings in turn lead to soft decay products, one can expect the high energies and boosts available at the LHC to be advantageous for probing many similar LLP models beyond the specific one discussed here. Such models might include, for example, inelastic-dark-matter models involving other coupling structures, such as the anapole or charge-radius operators [29][30][31].
In this connection, it is worth noting that other proposed experiments would also be able to probe parts of our region of interest within the parameter space of inelastic-dipole dark matter models. At CERN, these include the SHiP experiment [27,28], which would collide a 400 GeV proton beam with a fixed target. An estimate of the reach of this experiment within our parameter-space region of interest is provided in Appendix C. As can be seen there, the reach depends critically on the mass splitting between the dark states and the energy threshold for monophoton signal detection. For small ∆, FASER and FASER2 probe regions beyond the reach of SHiP, but for larger ∆, their reaches may be very roughly comparable. Further detailed studies of the energy threshold for monophoton signal detection are clearly needed. In Ref. [13], the authors explored the sensitivity of SHiP and FASER to an inelastic sterile neutrino interacting through an MDM operator and found SHiP to be more constraining than FASER. This is a consequence of their requiring very high thresholds for monophoton detection for FASER (E γ ≥ 10 GeV), while requiring much lower thresholds at SHiP (E γ ≥ 0.1 GeV). Although these thresholds are not precisely known for each experiment, FASER has significantly more shielding against background and is likely better suited to finding very lowenergy signal photons. We also note that other proposed experiments, such as SHADOWS [70] have the potential to probe this parameter-space region, although the same comments about energy thresholds apply for all fixed target experiments. Other LLP detectors such as ANU-BIS [71], CODEX-b [72], MAPP [73], MATHUSLA [74] are primarily designed to be sensitive to LLPs produced with significantly more transverse momenta than those considered here. Thus data from these detectors is not expected to be competitive with data from FASER in the parameter space we consider.
Finally, in this paper we have considered the case of a dark sector that comprises two states χ 0 and χ 1 that interact with each other and with the visible sector via an inelastic-dipole interaction. While this model is interesting in its own right, it can also be viewed as merely the simplest possible realization of a non-minimal dark sector consisting of N different states χ i with i = 0, 1, . . . , N −1 that interact with each other via a set of operators O ij , each with the MDM or EDM structure given in Eq. (2.1). Within this more general class of models, decay cascades involving multiple sequential LLP decays of the form χ i → χ j γ can develop. Scenarios in which decay cascades involving multiple LLPs can arise frequently give rise to distinctive phenomenological signatures [20] that differ from those observed in more traditional LLP scenarios. It would therefore be interesting to explore the prospects for probing dark-sector scenarios of this sort at FASER, FASER2, and other proposed experiments.
The squared matrix element for vector-meson decay to electrons is given by For the case of an MDM operator structure, the squared matrix element for the process V → χ 0 χ 1 is given by For the case of an MDM operator structure, the squared matrix element for this same process is given by We note that for m 1 = m 0 , the expressions in Eqs. (A3) and (A4) reduce to the corresponding results obtained in Ref. [31] for the case of a purely elastic dipole interaction. By contrast, the dominant contribution to χ 1 production from each relevant neutral pseudoscalar meson P arises is that associated with the three-body decay process P → γ * γ → χ 0 χ 1 γ, where γ * denotes an off-shell photon. The coupling between each P and a pair of photons is described by the effective Lagrangian where F P is the corresponding pseudoscalar decay constant and where F µν ≡ µνρσ F ρσ /2 is the dual electromagnetic field-strength tensor. The differential branching fraction with respect to the energy E * 1 of the final-state χ 1 particle in the rest frame of the decaying P can be written in the form where dΓ(P → χ 0 χ 1 γ)/dE * 1 denotes the corresponding differential partial decay width and where denotes the full partial width of P to photon pairs. This differential partial width is given by where |M(P → χ 0 χ 1 γ)| 2 is the squared matrix element for the decay process P → χ 0 χ 1 γ, summed over finalstate-particle spins. For the case of an MDM operator structure, this squared matrix element is given by where p, k γ , k 0 , and k 1 respectively denote the fourmomenta of P , γ, χ 0 , and χ 1 , and where q denotes the four-momentum of the virtual photon. By contrast, for the case of an EDM operator structure, this matrix element is given by We note that the factors of F P in Γ(P → γγ) and dΓ(P → χ 0 χ 1 γ)/dE * 1 cancel in Eq. (A6), while BR(P → γγ) is well approximated by its SM value Appendix B: Thermally-Averaged Cross-Sections and Decay Rates In this Appendix, we provide expressions for the relevant thermally-averaged quantities that appear in the collision terms Eq. (4.3).
In general, the thermally-averaged decay width for a particle that decays via the two-body process a → b + c, where the species c and d in the final state are in thermal equilibrium with the SM radiation bath, is given by where g a , n eq a , m a , and Γ a→bc , respectively denote the number of internal degrees of freedom, and equilibrium number density, mass, and proper decay width of the initial-state particle a. Since the equilibrium number density is where K ν (x) denotes the modified Bessel function of the second kind, one finds that The sole thermally-averaged decay width that appears in Eq. (4.3) is the one for the decay process χ 1 → χ 0 γ. The corresponding proper decay width from which this thermal average is derived is provided in Eq. (3.5).
The thermally-averaged cross-section for the 2 → 2 scattering process of the form a + b → c + d, where the species c and d in the final state are in thermal equilibrium with the SM radiation bath, can be written as [75] σv = 1 n eq a n eq b d 3 p a (2π) 3 where σ ab→cd is the invariant cross-section for the scattering process and where denotes the relative velocity of a and b in the background frame. To a very good approximation, when a and b are highly non-relativistic, this relative velocity coincides with the Møller velocity and the expression in Eq. (B4) simplifies to [76] σv = T 32π 4 n eq a n eq where s min ≡ (m a + m b ) 2 . For the case in which a and b are of the same species, this expression further reduces to We now turn to discuss the cross-sections for the individual processes whose thermal averages appear in Eq. (4.3). For both the MDM and EDM operator structures, the cross-section for the s-channel co-annihilation process χ 0 χ 1 → e + e − can be written as where s min = (m 0 + m 1 ) 2 and where the form of C (s) differs between the MDM and EDM cases. For the MDM case, this quantity takes the form where we have defined By contrast, for the EDM case, C (s) takes the form For processes in which χ 0 and χ 1 co-annihilate into hadronic final states, we model the inclusive cross-section using the R(s) function given in Ref. [77]. In particular, we take where Θ(x) denotes the Heaviside theta function and where m π is the pion mass. For both the MDM and EDM operator structures, the cross-section for the t-channel annihilation process χ 0 χ 0 → γγ can be written as where the form of C (t) differs between the MDM and EDM cases. In the MDM case, this quantity takes the form where we have defined in terms of the quantities By contrast, for the EDM case, C (t) takes the form where we have defined and where X ± are once again given in Eq. (B16). The corresponding expressions for the t-channel annihilation process χ 1 χ 1 → γγ in the MDM and EDM cases are identical in form to the expression in Eq. (B13), but with m 0 ↔ m 1 exchanged in this cross-section formula and all ancillary expressions. For both the MDM and EDM operator structures, the cross-section for the process χ 0 e ± → χ 1 e ± , in which a χ 0 particle up-scatters into a χ 1 particle off an electron or positron, can be written as σ 0e→1e = e 2 C (sc) The cross-section for the corresponding down-scattering process χ 1 e ± → χ 0 e ± may be obtained simply by exchanging m 0 ↔ m 1 everywhere in Eq. (B19) and all ancillary expressions. Moreover, the corresponding expression for the up-scattering of χ 0 into χ 1 off any other electrically charged SM particle may be obtained by replacing m e with the mass of that SM particle everywhere in Eq. (B19) and all ancillary expressions.

Appendix C: Prospects for Detection at SHiP
SHiP [27,28] is a proposed experiment at CERN in which a 400-GeV proton beam will collide with a fixed target consisting of tungsten and a titanium-zirconiumdoped molybdenum alloy. SHiP probes lifetimes comparable to FASER and FASER2 and so it is useful to compare the reach of these experiments.
As described in Ref. [78], SHiP is proposed to consist of a decay volume shaped as a 50-meter-long pyramidal frustum with an upstream face 45 meters from the target with dimensions 5 m × 2 m and a downstream face with dimensions 11 m × 6 m. The SHiP detector is expected to be able to detect photons with energies E γ 300 MeV [79]. However, resolving the monophoton signal that results from χ 1 → χ 0 γ decay from background may require a somewhat higher photon-energy threshold. For example, in assessing the prospects for detecting a monophoton signature of dark-photon decay in the context of a dark-axion-portal model, the authors of Ref. [80] adopted a threshold E γ > 2 GeV based on comparisons with the process in which the decay of the dark photon decay yields a e + e − pair rather than a photon. However, until a detailed monophoton study for SHiP is performed, it remains unclear precisely what threshold should be adopted in searches of this sort. Thus, in what follows, we consider a range of possible choices for this threshold.
The procedure we use for determining the expected number of signal events at SHiP for any given combination of m 0 , ∆, and Λ O is analogous to the procedure described in Sect. V B that we use in determining the expected number of events at other fixed-target experiments. We calculate the differential cross-section for π 0 production due to the scattering of protons in the beam with a molybdenum target using the BdNMC code package. We then scale this cross-section by N (POT) M as in Eq. (5.3) to obtain differential production cross-sections for heavier neutral mesons. The N (POT) M for all relevant meson species except the vector meson ψ(2S) were provided in Ref. [31]. We calculate the number of ψ(2S) mesons produced per proton using Pythia 8.2 and find that N (POT) ψ(2S) = 10 −6 . We then determine the expected number of χ 1 decays per proton on target within the decay volume of the SHiP detector from the branching fractions given in Appendix A, assuming the detector geometry described above.   We observe that SHiP has a comparable reach to FASER2 for large ∆, even for larger detection thresholds, but has far less sensitivity than FASER2 for small ∆, even for a detection threshold Eγ ≥ 0.3 GeV.
In Fig. 5, we show contours (orange dashed curves) of the discovery reach of the SHiP experiment within our parameter-space region of interest for 2 × 10 20 total protons on target -the number expected over roughly 10 years of operation [28]. Since this discovery reach depends sensitively on the threshold energy for a single photon, we plot contours for several different values of this threshold: E γ ≥ 0.3 GeV, 1 GeV, 3 GeV, and 5 GeV. As in Fig. 3, the panels shown in the left and right columns correspond to the choice of MDM and EDM operators, respectively, while the panels shown in different rows correspond to different values of ∆. Contours indicating the reach of FASER with an integrated luminosity of 300 fb −1 (red curves) and the reach of FASER2 with 3000 fb −1 (blue curves) are also shown for comparison.
We observe from Fig. 5 that for small values of ∆, the SHiP detector has little sensitivity within our parameterspace region of interest. This is primarily a reflection of the fact that for small mass splittings, the lab-frame energy of the photon produced by χ 1 decay typically lies below the threshold for detection -even if that threshold is as low as 300 MeV. By contrast, for large values of ∆, a larger fraction of the photons have energies that exceed the detection threshold and the discovery reach of SHiP is far greater. For example, for ∆ = 0.05, the reach of SHiP and FASER2 is very roughly similar for comparable running times.