Looking forward to inelastic DM with electromagnetic form factors at FASER and beam dump experiments

Inelastic Dark Matter (iDM) is an interesting thermal DM scenario that can pose challenges for conventional detection methods. However, recent studies demonstrated that iDM coupled to a photon by electric or magnetic dipole moments can be effectively constrained by intensity frontier experiments using the displaced single-photon decay signature. In this work, we show that by utilizing additional signatures for such models, the sensitivity reach can be increased towards the short-lived regime, $\gamma c\tau \sim O(1)\,$m, which can occur in the region of the parameter space relevant to successful thermal freeze-out. These processes are secondary iDM production taking place by upscattering in front of the decay vessel and electron scattering. Additionally, we consider dimension-6 scenarios of photon-coupled iDM - the anapole moment and the charge radius operator - where the leading decay of the heavier iDM state is $\chi_1 \to \chi_0 e^+ e^-$, resulting in a naturally long-lived $\chi_1$. We find that the decays of $\chi_1$ at FASER2, MATHUSLA, and SHiP will constrain these models more effectively than the scattering signature considered for the elastic coupling case, while secondary production yields similar constraints as the scattering.


I. INTRODUCTION
Dark Matter (DM) is an electrically neutral form of matter which, as a diverse set of observations indicates, makes up a significant portion of the energy content of the Universe [1][2][3][4][5].Arguably, one of the most motivated mechanism for DM production is its thermal freeze-out from the thermal plasma taking place in the early Universe [6,7].Null searches for Weakly Interacting Massive Particles (WIMPs) [8,9], which would acquire their relic density in this way, motivates the exploration of more elaborate scenarios.
An interesting proposal that shares many desirable features with WIMPs is the inelastic DM (iDM) [10].In its original form, it assumes that the Dark Sector (DS) consists of two states that are charged under a SM or dark U (1) gauge symmetry, and the coupling between the DS species and the gauge boson is non-diagonal.Moreover, although DM needs to be electrically neutral (it can be at most millicharged) [11,12], interactions with the photon can occur through electromagnetic (EM) form factors [13][14][15][16][17].
Combining these two ideas, the iDM scenario can be realized by introducing such a higher dimensional EM operator, e.g., an electric or magnetic dipole moment (EDM/MDM) [16,[18][19][20][21][22][23][24], while preserving the distinct inelastic character of the scattering signatures at collider or direct detection searches, as in the original scenario [10].One consequence of such coupling is the decay of the heavier iDM state into the lighter DS state and SM particles.If the masses of the iDM states are within ∼ sub-GeV range, and the mass splitting between them is small -incidentally, this regime is usually needed to obtain the correct relic density -the heavier state is a * k.jodlowski@ibs.re.kr long-lived particle (LLP).In fact, the iDM coupled to a dark photon is one of benchmarks for intensity frontier experiments looking for highly displaced LLP decays [25,26].
In particular, recent work has shown [27] that intensity frontier experiments such as the current FASER [28,29] and the upcoming FASER2 [30][31][32] and SHiP [33,34] detectors will cover much of the allowed parameter space for iDM connected to the SM via EDM or MDM.We note that the elastic coupling scenario was studied first, see, e.g., [35][36][37], but in this case the DM species is stable and therefore only its scattering signatures can be used.Moreover, since the LLP decays taking place inside a detector can be more efficient than scattering with electrons, it was found that the searches for decays of the heavier DS species, χ 1 → χ 0 γ, will cover substantially larger part of the parameter space than in the case of elastic coupling [27].
The presence of extended DS, that is two species with small mass splitting between them, allows one to search not only for LLP decays, but it also enables secondary LLP production [38][39][40].It takes place just in front of the decay vessel by upscattering of the stable DS species into the heavier one.
In this work, we extend [27] by considering secondary LLP production followed by its rapid decay and the electron scattering signature for EDM/MDM iDM.In addition, we consider additional scenarios for the iDM -when it is connected to the SM via anapole moment (AM) or charge radius (CR) operator.In these cases, the leading decay is phase-space suppressed, χ 1 → χ 0 e + e − , which could naturally explain the long-lived nature of χ 1 .We also study the impact of secondary LLP production and scattering with electrons.
The paper is organized as follows.In Sec.II, we introduce the scenarios for iDM connection to a photon via mass-dimension five and six operators.We review their properties, and in particular discuss the parametric de-pendence of LLP decay length, which can be probed at the intensity frontier.In Sec.III, we discuss the relic density of χ 0 obtained by thermal freeze-out.We identify the key annihilation/co-annihilation processes entering the Boltzmann equations, and we discuss their solutions.In Sec.IV, we discuss experimental signatures of long-lived iDM: displaced decays, secondary LLP production, and electron scattering in various experiments.We also provide details of our simulation.In Sec.V, we discuss our results, in particular the sensitivity plots of past and future experiments for the four iDM scenarios considered.We also compare our results with previous works.We conclude in Sec.VI, while expressions for the decay widths and cross sections used in our analysis were relegated to Appendices A to D.

II. MODELS
We consider iDM composed of two dark fermionsχ 0 , χ 1 -coupled to a photon by two dimension-5 operators, EDM and MDM, and by two dimension-6 operators, AM and CR operator.The key property of this scenario is that the DS states have split masses, ∆ ≡ (m χ1 − m χ0 )/m χ0 > 0.
This Lagrangian is an adaptation of the EFT framework used to describe fermionic DM with an elastic coupling to a photon [45,46] to iDM, keeping operators up to mass-dimension 6.We note that considering, e.g., scalar iDM is also possible, which we leave for further work.
Moreover, Eq. ( 1) is only an effective description valid for energies below the electroweak breaking scale, and in general requires UV completion, see, e.g., [20,21] for a discussion of such a construction for magnetic iDM and a dimension-7 Rayleigh operator.In fact, as dis-cussed in [16,18,23], Eq. ( 1) may arise from the elastic couplings case, where a Dirac fermion χ also obtains a Majorana mass, leading to a small mass splitting ∆.In such framework, another UV completion mechanisms have been discussed, e.g., [47][48][49][50][51]. Finally, replacing F µν with B µν , the hypercharge field strength tensor, substantially modifies phenomenology at, e.g., the LHC for m χ ≳ 100 GeV [52].Since we will be discussing a mass range of 1 MeV ≲ m χ0 ≲ 10 GeV, as well as the small momentum exchange regime in scattering processes involving the iDM species, Eq. ( 1) is an adequate description.
In this mass range, and for a small mass splitting, ∆ ≪ 1, the decay width of χ 1 is suppressed, and χ 1 acts as a LLP.Interestingly, the same regime is needed for successful thermal freeze-out of χ 0 from the thermal plasma in the early Universe thanks to the co-annihilation processes.This behavior is typical in iDM scenarios [10], and is discussed at length within our setup in Sec.III.
The lifetime of χ 1 strongly depends on the nature and dimensionality of the interactions.In particular, for EDM/MDM it is determined by two-body decays, χ 1 → χ 0 γ, with widths given by Eq. (A1).For dimension-6 operators, such two-body decays are not possible because the χ 0 -χ 1 -γ coupling depends on the four-momentum of the photon, which vanishes on-shell.On the other hand, three-body decays of χ 1 into χ 0 and a pair of charged leptons are allowed because the photon fourmomentum from the aforementioned coupling is canceled by the propagator of the off-shell photon.
Therefore, χ 1 → χ 0 e + e − , described by Eq. (A2), is the leading decay channel.Due to phase-space suppression, the χ 1 lifetime that is relevant for the intensity frontier searches in the AM/CR scenario is shifted towards larger values of m χ0 and ∆ compared to the EDM/MDM scenario.Moreover, we checked that this process does not lead to a noticeable injection of positrons relevant to the 511 keV excess in either decaying or excited DM scenarios (see discussion in the recent work [53] and references therein), since the χ 1 abundance (from freeze-out or upscattering χ 0 χ 0 → χ 1 χ 1 ) is negligible in the MeV-GeV mass regime we investigate.
where E is the energy of χ 1 the in the LAB frame, and we fixed all the parameters except the coupling constants to the typical values that can be covered by FASER2.We note that by using the above d χ1 formula for EDM/MDM case, we find disagreement with Eq. 2.3 of [27], which underestimates the lifetime by a factor of 260.Looking at the results shown in the bottom panels of Fig. 3 therein, the benchmark from Eq. 2.3 lies in the lower left part of the FASER2 sensitivity curve.However, this part of the sensitivity curve is determined by LLP species that are on the verge of being too long-lived to decay inside the detector, and therefore correspond to d χ1 ≫ 600 m.In the opposite regime, LLPs with a typical decay length d χ1 ≲ L, where L is the distance from their point of production to their decay (L ∼ 480 m for FASER) are in the upper right region of the sensitivity curve.Points above this line are too short-lived to be probed, since the probability of LLP decays taking place inside a distant detector is exponentially suppressed for d χ1 ≪ L [28,54].Therefore, the parameters used in Eq. 2.3 correspond to a much longer-lived LLP than indicated in [27].On the other hand, we believe that this discrepancy is a typographical error that does not impact any of their other findings.As discussed in Sec.V, we obtained similar sensitivity reaches for the EDM/MDM cases, which is described in that section's discussion.

III. RELIC DENSITY
Since the photon connects the visible and dark sectors, the stable iDM species χ 0 is a natural thermal DM candidate [10,16,21].The relic densities of AM/CR iDM2 is obtained by solving the following Boltzmann equations that describe the temperature evolution of the comoving yields Y i ≡ n i /s for i ∈ {χ 0 , χ 1 } ≡ {0, 1}: where λ ≡ π 45 We used: x = m 0 /T , r i = m i /m 0 , m Pl. = 1.9 × 10 19 GeV is the Planck mass, and g * s (T ) and g * (T ) are the effective number of degrees of freedom for the entropy and energy densities at temperature T of the SM-DS plasma, respectively.The brackets indicate the thermal average, e.g., ⟨Γ 1→0e + e − ⟩ = K1(m1/T ) K2(m1/T ) Γ 1→0e + e − , where K n (z) is the modified Bessel function of the second kind [55].
We numerically solve Eq. (3) using partial wave expansion [56] of the thermally averaged cross sections given by Eqs.(D1) to (D9).We obtain the present-day yields of each species i, Y i,0 , which are used to compute the resulting relic density where s 0 is the present-day entropy density.In all four iDM scenarios, the dominant process responsible for the successful thermal freeze-out is χ 0 -χ 1 co-annihilation into SM leptons and hadrons, ⟨σ 01→SM SM v⟩.Other processes, in particular χ 1 χ 1 → χ 0 χ 0 annihilation, are suppressed, this one by the fourth power of the coupling constant, and do not play a significant role in the thermal freeze-out.
As is typical for the co-annihilation mechanism [57], the correct relic density, Ω 0 h 2 ≃ 0.1, can only be obtained for sufficiently small values of mass splitting, ∆ ≲ O(0.5).This is precisely the regime of long-lived   3), while the dashed lines denote the equilibrium comoving yields.The orange dashed line indicates the relation given by Eq. ( 4), which is followed by Y1 over the shown range.
χ 1 , which allows to probe such iDM scenarios at the intensity frontier [23,27,58].This bound holds because in Eq. ( 3) the term corresponding to χ 0 -χ 1 coannihilations into leptons and hadrons is suppressed by a factor of ∼ exp(−x • ∆) compared to the same process with ∆ = 0.This suppression is due to the magnitude of the comoving yield Y 1 with respect to Y 0 ; note that the asymptotic behavior of the Bessel function at z → ∞ is K 2 (z) ∼ π/(2z) exp(−z) [55].
For the EDM/MDM case, we recover the results of [27].Therefore, in this discussion, we will focus solely on the AM/CR scenarios.In Fig. 1, we show the solutions of Eq. ( 3) for AM on the left and the solutions for CR on the right.We consider two mass splitting patterns, which will be discussed in more detail in Sec.V, corresponding to ∆ = 0.05 (top) and ∆ = 0.2 (bottom).
In both cases, iDM species decouple from the thermal plasma when the co-annihilation rate falls below the Hubble rate, which occurs near x ≈ 20.Note that the value of the coupling constant required to obtain the correct relic density is larger for the AM than for the CR operator due to the p-wave suppression of the co-annihilation cross section in the former scenario.
Moreover, similar to the case of EDM/MDM discussed in [27], once iDM departs from thermal equilibrium, the up-and down-scattering processes, as well as the χ 1 decays/inverse decays, ensure the validity of the following relation: This relation holds because in such co-annihilation freezeout scenario, chemical and kinetic equilibrium between iDM species is preserved even after the decoupling, thus then d(Y 0 + Y 1 )/dx ≈ 0, and since Y 1 ≪ Y 0 , one obtains Eq. ( 4).

IV. LONG-LIVED REGIME OF IDM
In this section, we describe the iDM production mechanisms and signatures probing the long-lived nature of χ 1 for the iDM scenarios introduced in Sec.II.
A. Intensity frontier searches for iDM a.Primary LLP production The production rates of iDM at FASER2 have been studied and discussed in previous works.For the elastic version of interactions given by Eq. ( 1), we refer to the discussion around Fig. 1 in [37], while for EDM/MDM iDM, we refer to [27] -Fig. 1 and the discussion there.
As discussed in these works (see also [36]), the main channel for the production of χ 0 -χ 1 pairs in the mass range m χ0 ≃ m χ1 ∼ MeV-GeV, which is relevant to the χ 1 long-lived regime, is vector meson decays, V → χ 0 χ 1 .Moreover, the branching ratio of such vector meson decays are roughly proportional to its mass squared [27,36].Therefore, the largest contribution comes from decays of the heaviest meson, provided there is sufficient abundance of it.The formulas for the branching ratio for all iDM operators are provided in Appendix B, see Eq. (B1).We also consider the next-leading contribution to the LLP production, which are three-body pseudoscalar meson decays, P → γχ 0 χ 1 .The double differential form of the branching ratio, convenient for Monte Carlo simulation, is given by Eq. (B2).
For both vector and pseudoscalar meson decays, we recovered results of [27] for EDM/MDM, and in the elastic limit (∆ → 0) of all four iDM operators, we found agreement with the results from [37], except for the pseudoscalar decay for MDM.In this case, we get the opposite sign in front of the sin 2 θ term (using their parameterization from Eq. A1 therein).We verified that in the m χ1 = m χ0 → 0 limit, when is proportional to the sin 2 θ term, the integrated formula from [37] yields a negative number.We note that this potential discrepancy does not affect any other results from [37] (or ours), since this decay is subdominant with respect to the vector meson decays.
b. Displaced LLP decays After the iDM species are produced, the heavier of the two, χ 1 , decays into χ 0 and a photon (a e + e − pair) via the EDM/MDM (AM/CR).The main signature considered in this work involves such decays taking place inside a dedicated detector, which is positioned at a significant distance (∼ 100 m) from the point of production of iDM species to mitigate the SM background.The corresponding probability of such an event is given by [59] p where d(E) represents χ 1 decay length in the laboratory reference frame, while L indicates the distance between the start of the detector, which has length ∆, and the LLP production point.As mentioned in the discussion following Eq.( 2), this probability is exponentially suppressed for LLPs characterized by insufficient decay length, d ≪ L, while in the opposite regime, d ≫ L, the suppression is only linear, p(E) ≃ ∆/d [60].Hence, L determines the minimal length scale for d that can be studied by such method.Intensity frontier searches looking for ∼ sub-GeV LLPs are typically only sensitive to decays to a pair of charged SM fermions or a pair of photons [25,26,61].Decays into a single high-energy photon are more challenging due to the significantly larger background compared to decays into e + e − or γγ.However, current LHC far-forward detector FASER [28,29,62] or beam dump experiments such as, e.g., the future SHiP [33,34] will have sensitivity to such signatures.In fact, several studies, e.g., [27,39,40,63], have already investigated various BSM scenarios with LLP decaying into a single photon.
c. Secondary LLP production As described above, searches using displaced LLP decays at remote detectors are limited to sufficiently long-lived BSM species.In addition to the main decay vessel of FASER, which is designed to study LLP decays, a specialized neutrino detector called FASERν [64,65] has been positioned in front of it.It is an emulsion detector made of tungsten layers, which allowed the FASER collaboration to detect collider neutrinos for the first time [66].
Moreover, FASERν can be used as a target to produce secondary LLPs in BSM scenarios with at least two DS species with different masses, which are produced in pairs in the forward direction of the detector. 3This process occurs through the upscattering of a lighter DS species into the heavier one, which acts as a LLP.It was shown in [38][39][40] that such upscattering is particularly efficient in the coherent regime, involving interactions with the entire nucleus (e.g., tungsten for FASERν).After the secondary LLP production, in our case χ 1 , it can decay within either the FASERν or FASER detectors. 4hus, in this way, FASER may cover part of the d ∼ L ≃ 1 m region of the parameter space that is inaccessible by primary production.Moreover, such a regime is particularly interesting, e.g., it typically corresponds  The upper panels correspond to E1 = 100 GeV of the incoming particle, χ0, while the lower panels correspond to E1 = 1000 GeV.Note the significant suppression of the cross sections for both processes for dimension-6 operators relative to dimension-5 operators due to the lack of 1/t amplitude enhancement.Since the cross sections given in Eq. ( 6) overlap for both EDM/MDM and AM/CR, we shifted the overlapping lines by 5%.For comparison, we also show the cross section for Primakoff upscattering for an axion-like particle (ALP) coupled to two photons.
to the correct thermal relic density region -see Fig. 6 in [38] for the case of iDM coupled to a dark photon, which connects to the SM via kinetic mixing.On the other hand, due to the form of the formula describing the secondary production cross section, see Eq. ( 6), the number of decays of LLPs produced in this way has an additional dependence on the coupling constant squared.Therefore, the secondary production of LLP can only be efficient for sufficiently large values of the coupling constant.
Secondary χ 1 production is dominated by coherent scattering with an entire nucleus due to the Z 2 enhancement, compared to only Z-dependent incoherent scatterings.We obtained the following integrated cross sections for such processes: where a = 111Z −1/3 /m e , d = 0.164 GeV 2 A −2/3 , m e is the electron mass, Z (A) is the atomic number (mass number) of a nucleus, and t max ≃ −(m 4 χ0 + m 4 χ1 )/(4E 2 1 ).We used the form factor given in Eq.B2 in [40].In addition, to obtain Eq. ( 6), we used the method described in [40,67].The full derivation, including the square of the amplitude, can be found in the accompanying Mathematica notebook.We consider two setups of the FASER2 experiment: the baseline scenario (left) and the Forward Physics Facility containing additional detectors, e.g., FLArE (right).Note the lower energy threshold (Eχ 0 γ > 1 GeV) for the photon coming from χ1 decay for the case ∆ = 10 −3 , which is needed due to kinematics of the process, to obtain a sufficient number of events.For both benchmarks, we reproduce results of [27] for primary LLP production (black and red solid lines for FASER and green dashed line for SHiP).Sensitivity reaches using secondary production (dashed and dash-dotted black lines) allow to cover the smaller lifetime regime, dχ 1 ∼ 1 m, which for ∆ = 0.05 largly overlaps with the thermal relic density line, Ωχ 0 h 2 = 0.1.The electron scattering signature at FASERν2 (solid gold line) and FLArE (dot-dashed gold line) leads to sensitivities analogous to the case of elastic coupling studied in [37].
d. Electron scattering signature In the case of elastic interactions described by Eq. ( 1), the DS species is stable, and the main signature at the intensity frontier experiments is elastic scattering with electrons [36,37].Following these works, we consider this signature for the iDM scenario.
The integrated cross sections are given in Eq. (C2), while Eq.(C1) gives their differential form.The latter is needed to impose the energy-and angle-dependent cuts for the FASERν2 and FLArE electron scattering searches [37,39,68]. 5Furthermore, following [39], we require that the χ 1 produced in the upscattering with electrons inside the FASERν decays outside of it and the main FASER decay vessel.In Fig. 3, we show the mass dependence of the cross section for upscattering (χ 0 T → χ 1 T , T = {N, e − }) on nucleus (left) and on electron (right) for each of the iDM operator we consider.We consider two energies of the incoming χ 0 , E 1 = 100 GeV and E 1 = 1000 GeV.We assumed a coupling constant for each model equal to 1 in units of GeV.Since the magnitude of the upscattering cross section plays a crucial role in the significant secondary production, we also present the cross section for the Primakoff production of ALP [67].For this model, the Primakoff process actually plays the role of primary LLP production, since it dominates the contributions of meson decays.However, for the extended ALP scenario with its coupling to a photon and a dark photon, called the dark axion portal [69,70], the analogous Primakoff process with the on-shell photon replaced by a dark photon has the same form, except it is smaller by a factor of 3/2 [40].As a result, such secondary production allows FASER2 to also cover part of the d ∼ 1 m region of the parameter space of this model, see Fig. 2-4 in [40].
In the iDM case, the character of the interactions is clearly visible -the dimension-5 operators lead to sizable cross sections, while results for operators of dimension-6 are suppressed.Therefore, one can expect to extend the sensitivity reach of FASER2 using displaced LLP decays thanks to secondary production only for dimension-5 operators.On the other hand, the cross sections given by Eq. ( 6) are enhanced by a factor of Z 2 , hence the number of such events may exceed the number of electron scattering, which is only proportional to Z. Thus, it is worthwhile to investigate the upscattering process acting as either secondary LLP production or scattering with electrons for all iDM scenarios.The origin of the suppression for AM/CR is similar to the reason why χ 1 → χ 0 e + e − decay is possible, while χ 1 → χ 0 γ is not. 6The χ 0 -χ 1 -γ coupling for the AM/CR operator depends on the four-momentum of the photon.Hence, the square of the amplitude vanishes for on-shell photon, while for off-shell photon, which mediates either the upscattering or decay of χ 1 , its propagator cancels this common factor of the coupling, rendering the square of the amplitude non-zero.This behavior also illustrates why the sensitivity reaches for electron scattering for elastic DS with EM form factors for EDM/MDM are much stronger than for the AM/CR operator [36,37].

B. Simulation of LLP signatures
We implement the signatures of the iDM operators introduced in previous subsection within modified FORESEE [71] package.The implementation includes, e.g., the production of LLPs as described in paragraph a) above.In particular, we used the SM meson spectra generated by EPOSLHC [72] and Pythia [73].We then use the routines of the modified FORESEE to obtain the χ 0 -χ 1 yield, and the number of events corresponding to each signature discussed in Sec.IV A.
We note that the experiments described above have been investigated for EDM/MDM iDM in [27], and we follow their discussion of simulating the LLP decays.In particular, we adapt their energy thresholds for all of these experiments.In addition, we examine other iDM signatures in the long-lived regime, which were described in the previous subsection.
Our choice of detectors is motivated by the desire to cover as much parameter space as possible.The selected experiments differ from each other in key physical characteristics, such as, e.g., the distance L, luminosity and energy of the proton beam, size, and visible energy threshold.As a result, as shown in the next section, they allow coverage of a variety of iDM scenarios, similar to many popular LLP scenarios, such as the renormalizable portals or an ALP coupled to a pair of photons; see extensive discussion in recent reviews [25,26,61].

V. RESULTS
In this section, we describe the prospects of detecting iDM with EM form factors in beam dump experiments and LHC far-forward detectors.All models described in Sec.II are characterized by three free parameters: (Λ, m χ0 , ∆) -the coupling constant, the mass of the stable iDM species, and the mass splitting.
As discussed in Sec.III, mass splittings ∆ plays a key role in determining the iDM relic density.Therefore, we present our results in the (m χ0 , Λ) plane for several values of ∆.In the case of EDM/MDM, we study the two benchmarks from [27], corresponding to ∆ = 10 −3 and ∆ = 0.05, respectively.On the other hand, for AM/CR, we investigate larger values of ∆: 0.05, 0.2, and 1.
This choice is motivated by the intention to examine the regime of long-lived χ 1 , which for dimension-6 operators, is shifted towards larger values of the mass splittings, cf.Eq. (A1) and Eq.(A2).We also note that for all iDM scenarios, our simulation is adapted to any ∆, and can be easily modified by, e.g., changing the experimental cuts.

A. Sensitivity reach for dimension-5 operators
In Figs. 4 and 5, we present the results for MDM and EDM iDM, respectively.For the displaced LLP decays produced by vector meson decays, which is the main iDM signature in the long-lived regime, we reproduce the results from [27].The legend of each plot includes the visible energy threshold imposed on the χ 1 decays in each experiment.The constraints shown were derived by LEP [84,85] and CHARM-II [86], while we also show the bounds from NuCal.In particular, note that the benchmark used in the first line of Eq. ( 2) is in the upper right area of the FASER2 sensitivity plot for both iDM operators, which agrees with our discussion following this equation.
For dimension-5 operators, our main results are the sensitivity lines obtained from the secondary LLP production.As described in Sec.IV A, this process proceeds through coherent upscattering, χ 0 N → χ 1 N , on the nuclei of the tungsten layers of the FASERν2 emulsion detector located in front of the main decay vessel, FASER2.Since the distance between the two detectors is ∼ O(1) m, secondary production allows to cover such short-lived LLP regime, which is otherwise challenging to probe because of the exponential suppression, see Eq. ( 5).We note that this configuration of the FASER detectors presents advantages over a typical beam dump experiment, in which there is a large separation, ∼ O(100) m, between the LLP production and detection points.
After χ 1 is produced by upscattering, it can decay in FASERν2 (dashed-dotted black lines) or inside FASER2 (dashed black lines).As can be clearly seen, this mode of χ 1 production allows covering the region of larger mass and coupling values than the primary production.This is particularly important for MDM, where it can help cover part of the parameter space where χ 0 is a thermal DM candidate.
We also present results for the electron scattering signature at FASERν2 (gold solid line) and FLArE (gold dot-dashed line).It covers the low-mass regime that complements the sensitivity reach obtained using decays of χ 1 produced in either primary or secondary production processes.The projection lines derived in such a way are generally weaker than those using secondary LLP production.This originates from the scaling behavior of the number of events with the atomic number, which is Z for the electron scattering and Z 2 for the upscattering production.We also note that the FASERν2 sensitivity ends at m χ0 ∼ 0.1 GeV due to the imposed condition on χ 1 decays -they must take place outside both FASER detectors -which is best illustrated in the bottom panels of Fig. 5.
In summary, secondary χ 1 production allows to cover a part of the short-lived regime of EDM/MDM that is inaccessible by primary LLP production.Depending on the value of the mass splitting ∆, it gives comparable or stronger bounds than the scattering signature, which for iDM gives similar results to the elastic DS with EM form factors scenario [36,37].As a result, a part of the allowed region of the MDM parameter space that corresponds to a successful thermal freeze-out of χ 0 will be covered by FASER2 during the high luminosity era of the LHC.

B. Sensitivity reach for dimension-6 operators
In Figs. 6 and 7, we show our results for AM and CR operator iDM.Both scenarios are mainly constrained by the bounds derived for the elastic coupling case, which were discussed in [36,37].These bounds were obtained using the results of LEP [84,85], SN1987A [87], and the LHC [52,88].We also show the exclusion plots obtained by us using NuCal.Moreover, similarly to the elastic case, in the AM/CR iDM scenario, the correct relic density of χ 0 is obtained for the parameter space that is already excluded.
As discussed in Sec.IV, the upscattering cross section for AM/CR is suppressed by several orders of magnitude compared to the EDM/MDM iDM, see Fig. 3.As a result, the secondary χ 1 production does not play a significant role, and the resulting sensitivity reach is comparable to the electron scattering signature, which in turn is similar to the scenario of elastic DS with EM form factors investigated in [37].
The main difference between that BSM scenario and ours is the possibility of χ 1 decays.As Eq. (A2) indicates, the leading decay width of χ 1 is phase-space suppressed, therefore, the lifetime regime that can be probed by intensity frontier searches is moved to larger values of ∆; see also discussion above Eq.( 2).On the other hand, such large mass splittings make the co-annihilation rate insufficient for successful freeze-out of χ 0 , leading to its overabundance.
To resolve this tension between the intensity frontier and cosmological considerations, in Figs. 6 and 7, we show the results for three values of the mass splitting ∆: 0.05, 0.2, and 1.In all cases, the displaced decays of χ 1 will probe a significantly larger parameter space than the electron scattering signature, whose reach is again similar to the scenario of elastic DS with EM form factors studied in [37].
Such behavior can be explained by the fact that the LLP decays are not suppressed by the electron number density like the electron scattering signature.Moreover, similar to the case discussed in Sec.V B, secondary production allows to cover a similar or larger area of parameter space than electron scattering.
On the other hand, for the largest value of ∆ we consider, the displaced χ 1 decays will allow FASER2, and in particular SHiP, to cover a part of the parameter space not covered by either the LHC or SN1987.Although in such a regime χ 1 cannot freeze-out to the observed DM relic abundance, it is still interesting that intensity frontier experiments will be able to extend the sensitivity reach beyond the LEP and LHC bounds.
In summary, the long-lived nature of χ 1 in the AM/CR scenario of iDM will allow LHC detectors such as FASER2 and MATHUSLA, and beam dumps like SHiP, to set significant exclusion limits that will be complementary to the existing bounds using other terrestrial searches.Moreover, in the regime of large mass splitting, such a signature will even allow going beyond such limits for a certain region of parameter space corresponding to m χ0 ∼ 1GeV.

VI. CONCLUSIONS
In this work, we have investigated the prospects of iDM with electromagnetic form factors at the FASER2, MATHUSLA, and SHiP detectors, while we also obtained bounds from previous experiments such as CHARM and NuCal.We focused on the long-lived regime of the heavier DS species, χ 1 , decaying into the lighter DS state and a photon (EDM/MDM portal) or a e + e − pair (AM/CR operator).
In the first case, we extended previous work studying χ 1 decays [27] by taking into the account the secondary LLP production occurring at FASERν2, which will allow FASER2 to cover the ∼ O(1) m region of the parameter space.This sector is particularly significant for MDM iDM, where in a part of this region, χ 1 is a thermal DM candidate.
For the AM/CR operator, we investigated for the first time the main iDM signature in intensity frontier searches, which is the displaced χ 1 decays.It allows FASER2 to be competitive with LEP bounds, with much greater reach than for the scenario of elastic DS with EM form factors, while SHiP may even cover a part of the allowed parameter space for large values of the mass splitting ∆.
In summary, we have shown that an iDM with electromagnetic form factors leads to greater detection prospects at FASER, MATHUSLA, and SHiP than the elastic coupling case, which was studied previously by DS-electron scattering signatures.The non-minimal iDM content allows not only for the decays of heavier DS species, but also for upscattering of the lighter one into the heavier one just in front of the decay vessel, which can lead to greater sensitivity than both displaced LLP decays and scattering with electrons.
Appendix A: Decays of χ1 In sections below, we collect the formulas used in our simulations.In some cases we provided only their leading form.The full results, including the derivations, can be found in the Mathematica notebook included in .
Below, we collect the decay widths for two-and threebody decays of χ 1 .For dimension-5 portals, we use results from [27]: (A1) It gives the expression for the photon energy E χ0γ , which is then used to impose the energy threshold cut.
For dimension-6 portals, such two-body decays are not possible, and the leading decay channel is the three-body decay of χ 1 into χ 0 and a pair of charged leptons.In the m χ0 ≫ m l limit, the width of this decay is (A3) where s 12 = (p 1 +p 2 ) 2 , s 23 = (p 2 +p 3 ) 2 , and p 1 , p 2 , p 3 are the momenta of χ 0 , l + , l − , respectively; the first factor is the double differential decay width, and the second one is the branching ratio of a virtual massive photon which accounts for the decays into hadrons taken from the PDG [89].
We also follow PDG kinematics chapter [89] for general expression for d 2 Γ ds12ds23 and the expressions for the integration limits.The amplitude squared averaged over the spins of the leptons reads: It is known that the two-body decays of vector mesons into χ 0 χ 1 pair dominate over the three-body pseudoscalar meson decays [27,36,37].Below we give the formulas for V (p 0 ) → γ * (p 1 +p 2 ) → χ 0 (p 1 )+χ 1 (p 2 ), which are mediated by an off-shell photon, EDM: where BR V→e + e − is the branching ratio corresponding to the V → e + e − decay, which we took from the PDG [89].

CR:
dBR P →γ χ0χ1 dq 2 d cos θ This form of the differential branching ratio allows for straightforward MC simulation of the pseudoscalar meson decays and is the form used in FORESEE.Below, we give the leading form of the expressions for the upscattering process, χ 0 T → χ 1 T , where the target is either an electron or a nucleus, T = {e − , N }.
where E R = E e − − m e − is the recoil energy of the electron.We have checked that in the ∆ → 0 limit we reproduce the elastic scattering results of [35].
The above expressions can be actually integrated analytically which results in EDM: Appendix D: Thermally averaged annihilation cross sections Below, we give all the the thermally averaged 2 → 2 cross sections entering into Eq.( 3).These are needed to determine the relic density of χ 0 , as described in Sec.III.We note that in ∆ → 0 limit we recover results of [35], when it is relevant.
Co-annihilation into charged leptons is given by Following [35], we take into account co-annihilations into hadrons according to the following formula: D3) where we use the data from PDG [89] for the experimentally measured R-ratio.
The annihilations into a photon pair is given by (D9)

Figure 1 :
Figure 1: Solutions of the Boltzmann equations shown for ∆ = 0.05 (top) and ∆ = 0.2 (bottom), with a fixed m0 = 0.5 GeV, resulting in the correct relic density Ω0h 2 ≃ 0.1.Results for the AM are on the left, while the results for the CR operator are on the right.The title of each plot contains the values of the parameters used to obtain each solution and the logarithm of the lifetime of χ1.The solid lines correspond to the numerical solutions of Eq. (3), while the dashed lines denote the equilibrium comoving yields.The orange dashed line indicates the relation given by Eq. (4), which is followed by Y1 over the shown range.

Figure 2 :
Figure 2: Contours of bχ corresponding to the observed DM relic density shown in the mχ 0 -∆ plane for the CR operator iDM.For anapole moment iDM, the resulting couplings are correspondingly larger, due to p-wave suppression.

Figure 3 :
Figure3: Mass dependence of the cross sections for upscattering taking place on tungsten (left) and electron scattering (right).The upper panels correspond to E1 = 100 GeV of the incoming particle, χ0, while the lower panels correspond to E1 = 1000 GeV.Note the significant suppression of the cross sections for both processes for dimension-6 operators relative to dimension-5 operators due to the lack of 1/t amplitude enhancement.Since the cross sections given in Eq. (6) overlap for both EDM/MDM and AM/CR, we shifted the overlapping lines by 5%.For comparison, we also show the cross section for Primakoff upscattering for an axion-like particle (ALP) coupled to two photons.

Figure 4 :
Figure4: Sensitivity of FASER2 to inelastic MDM for ∆ = 10 −3 (top) and ∆ = 0.05 (bottom).We consider two setups of the FASER2 experiment: the baseline scenario (left) and the Forward Physics Facility containing additional detectors, e.g., FLArE (right).Note the lower energy threshold (Eχ 0 γ > 1 GeV) for the photon coming from χ1 decay for the case ∆ = 10 −3 , which is needed due to kinematics of the process, to obtain a sufficient number of events.For both benchmarks, we reproduce results of[27] for primary LLP production (black and red solid lines for FASER and green dashed line for SHiP).Sensitivity reaches using secondary production (dashed and dash-dotted black lines) allow to cover the smaller lifetime regime, dχ 1 ∼ 1 m, which for ∆ = 0.05 largly overlaps with the thermal relic density line, Ωχ 0 h 2 = 0.1.The electron scattering signature at FASERν2 (solid gold line) and FLArE (dot-dashed gold line) leads to sensitivities analogous to the case of elastic coupling studied in[37].
Appendix C: Cross sections for electron orPrimakoff upscattering
for ⟨σv⟩ χ1χ1→γγ is obtained from the above by replacing each m χ0 with m χ1 and vice versa.The χ 1 -χ 0 conversion processes are determined by pairEM m e ∆m χ0 m e , CR:C = 3 √ 2α EM b 2χ m e ∆m χ0 m e .
23(m χ1 +m χ0 ) 2 −2s 12 +s 2 23 .In the case of three-body, we impose the energy threshold cut in the following way.For each point in the (Λ, m χ0 , ∆) parameter space, we calculate the average visibleenergyof the e + e − pair as a function of E χ1 , (A5) ⟨E e + e − ⟩(E χ1 ) ≡ ds 12 ds 23 d 2 Γ χ1→χ0e + e − ds 12 ds 23 × E e + + E e − Γ χ1→χ0e + e − , where we use PDG [89] for expressions for energies of e ± , E e ± ≡ E e ± (E χ1 , s 12 , s 23 ).Next, for a given energy threshold, E th , we determine E χ1 for which ⟨E e + e − ⟩ > E th .Finally, we use this value of E χ1 as a cut on χ 1 .