Searching for inelastic dark matter with future LHC experiments

We consider a dark sector containing a pair of almost degenerate states coupled to the Standard Model through a dark photon mediator. This set-up constitutes a simple realization of the inelastic dark matter scenario. The heaviest dark state is long-lived, in the limit of a small kinetic mixing among the dark photon and the Standard Model hypercharge gauge boson, and/or of a small mass splitting among the dark states. We study the prospects for detection of this scenario at proposed LHC experiments dedicated to search for long-lived particles, namely FASER, MATHUSLA, CODEX-b, AL3X, MAPP, ANUBIS and FACET. We consider both the cases of fermionic and scalar inelastic dark matter. We show that these experimental facilities can probe unexplored regions of the parameter space of this model, and we highlight their complementary roles.


Introduction
In recent years, a rich experimental program has been put forward to search for Long Lived Particles (LLPs) [1,2]. In accelerator experiments LLPs travel macroscopic distances before decaying, leading to signatures such as displaced vertices or missing energy. To fully exploit the potential of the Large Hadron Collider (LHC), new experimental facilities dedicated to search for LLPs have been recently proposed, namely FASER [3,4], MATH-USLA [5][6][7], CODEX-b [8,9], AL3X, [10], MAPP [11,12], ANUBIS [13] and FACET [14]. These proposed detectors would be located near the LHC Interaction Points (IP), and are designed to detect the decays of LLPs produced at the LHC, having at the same time the capability to dramatically reduce backgrounds. Interestingly, these projects are highly complementary to the existing LHC experiments ATLAS, CMS, and LHCb.
LLPs are predicted in many extensions of the Standard Model (SM) of particle physics, see e.g. [6,15,16]. In particular, they often arise in dark sectors, which are scenarios postulating the existence of new states, possibly including Dark Matter candidates (DM), feebly coupled to the SM. In this work, we focus on the possibility of new SM singlets in the form of inelastic DM (iDM) [17]. This model contains two particles almost degenerate in mass, the lightest one being the DM candidate (of course in general an extension with more than two states is possible). If this pair of particles is then coupled to the SM only through tiny portal interactions, the decay of the heaviest dark state is suppressed and this state acts as a LLP. iDM is particularly attractive from the phenomenological point of view, because it allows to explain the cosmological DM abundance via the standard thermal freeze-out and, at the same time, evade existing bounds. Coannihilations of the iDM particles in the early Universe determine their relic abundance. On the other hand, these processes are suppressed at later times, allowing to fulfill the constraints from indirect DM searches and Cosmic Microwave Background (CMB) observations, which are particularly tight for DM masses few O(10) GeV. The inelastic scattering of DM into nuclei is also suppressed, for large enough mass splitting among the dark states. Therefore, also direct detection constraints can be accommodated in this scenario. In general, the portal interactions among the dark sector and the SM can be described in terms of renormalizable (see [6,15,16] and references therein) or non-renormalizable portal operators [18][19][20]. We adopt the first option, considering a mediator between the SM and the dark sector in form of a dark photon [21,22]. A relatively light mediator is required to explain the DM abundance via standard thermal freeze-out processes. Building upon ref. [21], we perform an analysis of this iDM model, focusing on future searches of LLPs at the LHC. Our goal is to extend the analysis of [21] to recently proposed experimental facilities 2 . In particular, we determine the region of parameter space that can be explored by these proposed future LHC detectors, highlighting their complementary role in covering the parameter space of iDM models.
The paper is organized as follows. In Section 2 we introduce the iDM models under consideration. We will consider both scalar and fermionic iDM states. In Section 3 we describe how the projected sensitivities are computed. We present our main results in Section 4 and we offer our conclusions in Section 5.

Inelastic dark matter
Models of iDM have been proposed long ago as a simple way to evade the strong constraints coming from direct detection experiments [23,24], and have later been used to reconcile the observations by the DAMA experiment with the exclusion limits derived from other direct detection experiments [17] 3 . Although recent analysis seems to exclude an iDM interpretation of the DAMA signal [25], in recent years iDM has been studied as a paradigmatic example of a dark sector that can evade the stringent bounds coming from direct detection [26] and CMB observations [27], with DM particles still being thermally produced.
To illustrate the main idea, let us consider a scenario where, at leading order, the interaction between the dark and SM sectors occours via a trilinear vertex of the form D 1 -D 2 -mediator, where D 2 is a dark state obeying m 2 m 1 (with m 1 and m 2 the D 1 and D 2 masses, respectively) and D 1 is the DM candidate. The mediator involved in this vertex is coupled to SM particles, providing therefore a portal between the dark sector and the SM. In an underground laboratory experiment, for a sufficiently large mass splitting the incoming DM particle D 1 does not have enough energy to upscatter off a nucleus into the heavier state D 2 . Therefore bounds from direct detection experiments are evaded. At the same time, the iDM model can also fulfill the constraints on DM annihilations derived from CMB data. In fact, in the early Universe, the cross section for co-annihilation of iDM pairs D 1 D 2 into SM particles is suppressed by an exponential term exp[−(m 2 − m 1 )/T ] [28]. These co-annihilation processes are therefore irrelevant at late times, after the recombination era, as long as (m 2 − m 1 )/T rec 1, where T rec is the temperature at which recombination happens. For weak-size interactions, any mass splitting O(10) eV is sufficient to avoid the CMB bounds. More generically, for sufficiently large mass splitting, the constraints from indirect detection of DM are fulfilled in iDM models. On the other hand, as long as the mass splitting is not too large, the exponential factor is of O(1) at temperatures T ∼ m 1 , allowing for DM thermal production in the early Universe.
Relevant constraints might arise from annihilations of D 1 pairs, for instance processes such as D 1 D 1 ↔ ff , where f are SM fermions. Similarly, also interactions mediating the D 1 elastic scattering off nuclei might be present. These contributions are model dependent and will be discussed below.
As we will see in a moment, the situation presented above can be realized in a simple way for both fermionic and scalar dark states. For definiteness, we will consider the case in which the mediator between the dark sector and the SM is the so-called dark photon [29], i.e. the (massive) gauge boson associated with an additional U (1) symmetry, that interacts with the SM via a kinetic mixing term: In the previous equation, B µν and A µν are the field strengths of the gauge bosons associated with the hypercharge U (1) Y and the additional U (1) , respectively, cos θ w is the cosine of the weak angle and is the kinetic mixing parameter. In the remainder of the paper, we will remain agnostic about the origin of the dark photon mass, supposing that any scalar degree of freedom associated with the spontaneous breaking of U (1) is sufficiently decoupled from the D 1 , D 2 and A system. In the limit 1 and m A m Z (with m A and m Z the dark photon and Z boson masses, respectively), all SM fermions acquire a coupling to the dark photon of the form [30] where e Q f is the electric charge of the SM fermion f . This corresponds to the limit of photon-like interactions of the dark photon. However, since in our analysis we scan over the dark photon mass, considering also cases with a significant mixing between the A and Z bosons, we use the full expressions for the couplings of A to the SM states and of Z to the dark states [30].

Fermionic case
We now turn to a brief discussion of two realizations of the inelastic dark matter scenario described above. We start with the case of fermionic dark states. Let us take two Weyl fermions ψ 1 and ψ 2 with charges +1 and −1 under U (1) , respectively. The relevant Lagrangian besides standard kinetic terms (and written using two-component spinors) is given by where g d is the gauge coupling corresponding to U (1) and the Majorana masses δm 1,2 explicitly break U (1) . The latter can be generated, for instance, by the vacuum expectation value of some scalar state (in the minimal model, the same one that generates the dark photon mass). Notice that, since in the δm 1,2 → 0 limit the U (1) symmetry is recovered in the fermionic sector, it is technically natural to have δm 1,2 m D . In this limit the heavier state will typically be long lived. Taking all the mass parameters real for simplicity, the mass eigenvalues at leading order in corresponding to the mass eigenstates The factor of i in the definition of χ − amounts to a choice of phase that automatically guarantees m 1 > 0 under the assumption m D > 0. With this choice of phase, the interactions with the dark photon can be expressed in terms of Majorana spinors χ 1 = (χ − , χ † − ) and χ 2 = (χ + , χ † + ), obtaining This is precisely a vertex of the form D 1 -D 2 -mediator described above, with D 1 = χ 1 and D 2 = χ 2 . The terms suppressed by (δm 1 − δm 2 )/m D contain diagonal couplings of the typeχ i γ µ γ 5 χ i A µ , with i = 1, 2. In the limit of photon-like couplings as in eq. (3), the annihilation cross-section for the process χ 1 χ 1 ↔ ff induced by these terms is of p-wave type, and thus it is typically sufficiently suppressed to fulfil the CMB bounds mentioned above. In the following, we will assume that these terms are negligible or simply not present (implicitly taking the limit (δm 1 − δm 2 ) → 0). Notice that loop processes can induce χ 1 elastic scattering off nucleons, as well as χ 1 pair annihilations. These processes are however suppressed, below the sensitivity of current probes [21].

Scalar case
A similar reasoning can be applied to the case in which the dark sector is populated by scalar states. Considering a complex scalar φ = (φ 2 + iφ 1 )/ √ 2 with charge +1 under U (1) , the relevant Lagrangian contains the terms: where the mass term δm 2 explicitly breaks U (1) and splits the real scalars φ 1 and φ 2 , which acquire squared masses m 2 ∓ δm 2 , respectively. As in the previous case, it is technically natural to have δm 2 m 2 and for φ 2 to be long lived. As for the interactions with the dark photon, in terms of the mass eigenstates φ 1,2 they become once more of the form D 1 -D 2 -mediator with D 1 = φ 1 and D 2 = φ 2 . In this case, the annihilation cross-section φ 1 φ 1 ↔ ff , as well as the φ 1 elastic scattering off nuclei, are generated at the 1-loop level and will thus be parametrically suppressed by a factor ( /(4π)) 4 , sufficient to avoid the present constraints.
In what follows we will use the interaction Lagrangians of Eqs. (7) and (9) for our phenomenological analysis. In the discussion above and in the following, we consider the case of a dark photon heavier than the iDM pair, i.e. m A > m 1 + m 2 . The opposite limit corresponds to the secluded dark matter scenario [31], which has a different phenomenology. In that case, at colliders, which is the main focus of our work, A will decay into SM particles, rather than iDM pairs. Moreover, in the early Universe, the DM abundance is controlled by the annihilations D 1 D 1 → A A . These processes are strongly constrained by CMB and indirect detection bounds, see e.g. [32,33].
Before concluding this section, we discuss current experimental bounds on the dark photon parameter space. Depending on the relative importance of and g d , the dark photon may decay mainly visibly (into SM particles) or invisibly (into dark sector particles) [34]. Since we are interested in the case in which the interactions between A and the dark pair D 1 D 2 are sufficiently strong to allow for DM thermal production, in this work we consider only the so-called "invisible" dark photon. Focussing on the region m A 1 MeV, the most stringent bounds come from the dedicated search at BaBar [35] and electroweak precision measurements at LEP [30]. The BaBar limits dominate in the region m A 8 GeV, constraining 10 −3 in most of the parameter space and becoming as strong as 4 × 10 −4 for masses m A ∼ (4 ÷ 6) GeV. On the other hand, LEP bounds dominate for larger A masses and, as long as m A = m Z , put a milder bound of 2 × 10 −2 . Around m A m Z the limit becomes stronger, improving by roughly one order of magnitude. In addition, in this region the mixing between the A and Z bosons is close to maximal, with the Z boson potentially decaying into the dark states. For this reason, we also consider the bounds from Z → invisible [36].

Experimental setup and sensitivity estimation
Our analysis focuses on a scenario with (fermionic or scalar) iDM and with a dark photon mediator heavier than the pair of dark states, i.e. m A > m 1 + m 2 . As explained in the previous section, this construction allows the evasion of the stringent bounds from direct detection, indirect DM searches, and CMB observations. Furthermore, we are going to consider masses of the dark sector particles heavier than few GeVs. As we will see later, this region is favored by considerations on the DM relic abundance. The dominant production mechanism of iDM at LHC are Drell-Yan processes with a A or Z boson in the s-channel. An additional channel for the production of A particles, which would then decay into iDM pairs, is bremsstrahlung. This is relevant for m A few GeVs, and it is therefore subdominant in the range of masses than we explore. We implement the iDM models in Feynrules [37] modifying the Feynrules model file of [30], and use MadGraph5_aMC@NLO [38]

Forecasting sensitivity
The number of signal events in the detectors introduced in Sec. 1, can then be computed by multiplying the total number of D 2 particles produced at the LHC, by the probability that the D 2 particles decay inside the detector and lead to signal events, f signal . The first term can be simply obtained from the production cross-section of D 1 D 2 pairs, that we evaluate with MadGraph5_aMC@NLO, and the total integrated luminosity received by a given experiment, listed in Sec. 3.2. Instead, the quantity f signal can be computed as: where f dec corresponds to the probability for a D 2 particle to decay inside the detector and it is given by: It depends on the distance between the LHC IP and the point at which the D 2 particle enters (exits) the detector L entry (L exit ), and the decay length of D 2 in the laboratory (LAB) frame, given by Here τ D 2 is the decay time of D 2 , while β D 2 and γ D 2 are respectively its speed in units of speed of light (c), and its Lorentz factor, in the LAB frame. Clearly, f dec vanishes for those trajectories of D 2 which do not intersect the detector. In the region of interest, the main contribution to the D 2 decay width comes from the decays D 2 → D 1 ff into a pair of SM fermions. While the contribution into leptons is straightforward to compute, for quarks we need to worry about the validity of perturbative QCD. We follow the strategy outlined in [19] and use a perturbative computation when m 2 − m 1 > 2 GeV, while in the opposite case we use chiral perturbation theory and vector mesons matrix elements to capture the hadronic decays. For more details, see [19]. Eq. (10) is computed averaging ( · ) over all the possible kinematical configurations of the D 2 particles. For this purpose, we produced a sample of D 1 D 2 events with MadGraph5_aMC@NLO, that we then use to statistically evaluate eq. (10). Finally, in eq. (10), the quantity det takes into account the efficiency for the reconstruction of the events (that we simply take as 100%), and selection cuts. For the latter we impose a requirement on the energy of the visible final states employing the following strategy. Using MadGraph5_aMC@NLO we generate a sample of decays of D 2 → D 1 ll, where l is a lepton. We then perform the appropriate Lorentz boost to go from the D 2 frame to the LAB frame, we evaluate the sum of energies of the two leptons E vis , the visible particles, and we impose E vis > E cut . The fraction of the events satisfying this requirement leads to the quantity det in eq. (10). This procedure can be used to evaluate det for each D 1 D 2 event produced in our simulation. For FASER and FASER 2, we take E cut = 100 GeV, mimicking the cut considered in [3], for AL3X we adopt E cut = 3 GeV, as in [10], while for FACET E cut = 10 GeV [14]. For all the other experiments, we take E cut = 1 GeV, which is in the range of energy thresholds considered in [6]. In the right panel of Fig. 2 we show how the sensitivity changes doubling these reference values of E cut , for the representative cases of FASER 2 and MATHUSLA, respectively forward and off-axis detectors. We compute projected 95% C.L. limits imposing that the number of signal events is larger than N = 3. With this procedure we are assuming that backgrounds in the different experiments can be reduced at a negligible level, see e.g. the discussions in [3,4,6]. Finally, we compute the DM relic abundance using the public tool micrOMEGAs5.2 [39].

Experimental geometries
We shall now describe the experimental facilities that we consider in our analysis. The geometry of the MATHUSLA experiments is detailed in [7]. The detector is delimited by: where the coordinate system is centered at the LHC interaction point, the z axis is along the beam direction, and x denotes the vertical to to the surface. For a given trajectory of the χ 2 particle from our simulation, we evaluate whether the trajectory intersects the detector, and computes the quantities L entry and L exit . The total expected integrated luminosity is 3 ab −1 . We then consider the forward FASER detector and its proposed upgrade FASER 2 [3,4]. The FASER detector will be a cylinder with a radius of R = 10 cm, a length of 1.5 m, located 480 m from the ATLAS IP. The relevant integrated luminosity will be 150 fb −1 . For FASER 2 we take R = 1 m, a length of 10 m, a distance of 650 m from the IP, and an integrated luminosity of 3 ab −1 , [40,41].
Concerning ANUBIS [13], a cylindrical detector displaced from the interaction point is proposed. For simplicity, we approximate the circular base of the cylinder with a square of equal area. We take an integrated luminosity of 3 ab −1 . Another implementation of the ANUBIS geometry is provided in [42], see also [43]. We checked that using that parametrization leads to sensitivity curves consistent with our results.
The proposed CODEX-b would be installed next to the LHCb IP [8,9]. Its geometry is defined by and the total luminosity would be 300 fb −1 .
We then consider the AL3X detector, proposed to be placed close to the ALICE experiment at the LHC [10]. It is a 12 m long cylinder centered around the beam axis about 4.25 m from the IP with inner and outer radii respectively of 0.85 m and 5 m. This forward detector will cover the pseudorapidity range 0.9 η 3.7. We take an integrated luminosity of 250 fb −1 , which correspond to the optimistic value proposed in [10].
The MAPP-1 and its upgrade MAPP-2 detectors will be placed in the UGCI gallery adjacent to the MoEDAL region [11,12]. The projected integrated luminosity for these facilities is respectively of 30 fb −1 and 300 fb −1 . Concerning their geometry, we adopt the same setup of [42].
Very recently, a forward detector, FACET, to be placed around the CMS interaction point, has been proposed [14]. The detector volume is 18 m long, covering 101 m < z < 119 m, and it is sensitive to particles with polar angles 1 mrad < θ < 4 mrad. The expected total luminosity received is 3 ab −1 .

Results
In this section we present the results of the forecasted sensitivities using the methods and experiments detailed in Sec. 3. These are shown in the m 1 − plane for DM masses m 1 ≥ 1 GeV. This region is interesting because ample areas are unconstrained by current experiments and, as we are going to see, the LHC experiments listed in sec. 3.2 can probe large portions of it. In the same region we can also have thermal DM production. However, we do not restrict only to the parameters that allow to obtain the measured DM thermal abundance, since we can imagine a non-thermal DM production or a modified cosmological history that could dramatically change the picture. In all our plots we fix α D = g 2 d /(4π) = 0.1. To illustrate the complementarity of the different future LHC experiments, we focus on four representative benchmark scenarios, considering two different values for the mass splitting (∆ = 0.1 and 0.01) and for the ratio between the dark photon and the DM mass (m A /m 1 = 3 and 6). The main results are presented in Fig. 1 and Fig. 2. The shaded grey region depicts the current experimental constraints on the invisible dark photon already mentioned in Sec. 2, coming from BaBar [35] and LEP [30]. Since in the region m A ∼ m Z the Z boson coupling to the dark states is not suppressed, we also include the bound coming from the Z invisible decay [36]. The colored contours, on the other hand, show the projected future sensitivities of the experiments listed in Sec. 3.2. Finally, along the dashed black line the D 1 abundance matches the observed one. Fig. 1 and Fig. 2 correspond to the fermionic iDM (see Sec. 2.1) and scalar iDM (see Sec. 2.2) models, respectively.
The forecasted sensitivities are very similar between the fermionic and scalar cases for most of the experiments. This can be understood in the following way: on the one hand, the lifetimes of χ 2 and φ 2 are very similar, and the same is true for production of fermionic and scalar iDM pairs from on-shell dark photon decays. On the other hand, events with very off-shell dark photons give rise to different energy distributions depending on the dark states spin, and this explain the differences observed among the two scenarios. For this reason, in Fig. 2 we only show the benchmark ∆ = 0.1 and m A /m 1 = 3. Concerning the predictions for the DM relic abundance, the difference between the fermionic and scalar cases is instead more pronounced, due to the different velocity dependence of the non-relativistic annihilation cross-section.
As shown in Fig. 1 and Fig. 2, we find that most of the experiments will be able to probe regions of parameter space not excluded by current data and in which D 1 can be 100% of the DM. More specifically we see that, for ∆ = 0.1, the thermal line will be completely probed for m 1 10 GeV, while for ∆ = 0.01 the thermal line will be approximately covered in all the parameter space considered.
Comparing the ∆ = 0.1 and 0.01 cases, we see that for the latter choice the sensitivity of LHC experiments moves to larger values of m 1 and . This can be qualitatively understood looking at the decay width of D 2 → D 1 ff . In the limit of massless fermions and for small ∆ we have Γ ∝ 2 ∆ 5 m 5 1 /m 4 A [21]. A decrease of ∆ can be compensated with an increase of and m 1 to obtain the same decay width (and hence lifetime), reproducing the behaviour observed. Instead, moving from ∆ = 0.1 and 0.01 has only a modest impact on the predictions for the DM relic abundance. The situation is different for what concerns the role of the ratio between the dark photon and the DM mass. Going from m A /m 1 = 3 to m A /m 1 = 6 reduces the impact of the dark photon resonance for the DM annihilations in the early Universe. This implies that larger values of should be considered to obtain the same DM relic abundance. Let us also stress once more that the predictions for the DM abundance could be completely different in presence of non-thermal production mechanisms or scenarios of non-standard cosmology.
Concerning the different experiments, as evident from the discussion in Sec. 3, the sensitivity depends on a combination of several ingredients: the proximity of the detectors to the IP, their off-axis vs on-axis geometry, their size, luminosity received, and other factors such as the selection cuts performed in the analysis. As shown in Fig. 1 and Fig. 2, the different facilities considered here are complementary, but at the same time they cover overlapping regions of the parameter space.
The sensitivity reach of the AL3X experiment is particularly remarkable, although the feasibility of this proposal strongly depends on the overlap with the ALICE physics program [10]. The ANUBIS experiment allows to reach sensitivities comparable to the ones for MATHUSLA despite the smaller detector size, a factor which is compensated by its proximity to the IP. From Fig. 1 and Fig. 2 one can notice that MATHUSLA typically probes slightly smaller values of , but ANUBIS allows to test complementary parts of the parameter space (see e.g. the region around m 1 ∼ 20 GeV in Fig. 2). The recently proposed forward detector FACET also offers good sensitivities. The main advantage with respect to the FASER 2 is its proximity to the IP (101 m for FACET and 650 m for FASER 2) and the larger solid angle coverage. These features allow an increased sensitivity with respect to FASER 2, at least as long as the background can be reduced to negligible levels as assumed in the original proposal [14]. The projected reach of MAPP 2 is complementary to the ones previously mentioned, despite the smaller integrated luminosity that this detector will receive (0.3 ab −1 instead of 3 ab −1 ). At masses m 1 10 GeV, the sensitivity curves are qualitatively similar to the ones for MATHUSLA, but shifted to larger values of . Overall, these facilities are complementary to CODEXb, FASER (FASER 2) and MATHUSLA, previously considered in [21], surpassing their sensitivities in some regions of the parameter space. At qualitative level, similar results have been found for different models, see e.g. [14,[42][43][44]. However, the comparison of one particular experiment over the others strongly depends on the physical scenario under consideration, motivating therefore our dedicated study for iDM. Concerning CODEXb, FASER 2 and MATHUSLA, we have found results similar to [21], once the same setup is adopted. For this study, we have used the most updated geometry proposals for FASER 2 and MATHUSLA. These lead to slightly decreased performances for the former detector, and increased sensitivities for the latter (especially around m 1 20 GeV for the scenarios with ∆ = 0.1 and m A /m 1 = 3).
Finally, we shall mention that searches of iDM complementary to those analyzed here can be performed at ATLAS, CMS and LHCb looking for displaced muons and time-delayed tracks [21], and at Belle II with displaced vertices and missing momentum signatures [45,46].

Conclusions
We have considered a model of fermionic (and scalar) iDM coupled to the SM via a dark photon mediator. In this simple scenario, even for relatively light DM, it is possible to obtain the observed cosmological DM abundance via the standard freeze-out mechanism and, at the same time, fulfill the constraints from direct and indirect DM searches and CMB observations. The heaviest iDM state can be long-lived and it can be searched for at proposed LHC experiments. In this paper, we have updated and extended the analysis in [21], considering the future detectors FASER, MATHUSLA, CODEX-b, AL3X, MAPP, ANUBIS and FACET. The main results are shown in Fig. 1 and Fig. 2. We find that the experimental facilities discussed here offer promising prospects for detection. They can cover complementary regions of the parameter space, significantly extending the reach of current and past experiments. In the case of the mass ratio m A /m 1 = 6, and for the range of ∆ considered 0.01 < ∆ < 0.1, they will be able to completely probe the parameter space in which the lightest state constitutes 100% of the DM via the standard thermal freeze-out scenario.