Dark matter interacting via a massive spin-2 mediator in warped extra-dimensions

We study dark matter interacting via a massive spin-2 mediator. To have a consistent effective theory for the spin-2 particle, we work in a warped extra-dimensional model such that the mediator(s) are the Kaluza-Klein (KK) modes of the 5D graviton. We pay close attention to dark matter annihilations into KK-gravitons. Due to the high energy behavior of longitudinal modes of spin-2 fields, these channels exhibit a tremendous growth at large center of mass energies $\sqrt{s}$ if only one spin-2 mediator is considered. For the first time, we include the full KK-tower in this dark matter production process and find that this growth is unphysical and cancels once the full field content of the extra-dimensional theory is taken into account. Interestingly, this implies that it is not possible to approximate the results obtained in the full theory with a reduced set of effective interactions once $\sqrt{s}$ is greater than the first graviton mass. This casts some doubt on the universal applicability of previous studies with spin-2 mediators within an EFT framework and prompts us to revisit the phenomenological allowed parameter space of gravitationally interacting scalar dark matter in warped extra-dimensions.


Introduction
Over the last decades, a stupendous amount of evidence for the presence of Dark Matter (DM) in our Universe has been collected. This observational program culminates in the precision measurement of the CMB anisotropies by the Planck collaboration that allows a percent level determination of its abundance [1]. However, all these observations only test the gravitational effect of DM and despite large efforts, a detection of other interactions of the DM with the Standard Model (SM) is still elusive. As interactions mediated by SM particles and popular SM extensions are increasingly constrained, see e.g. [2][3][4][5], alternative ideas are hotly sought after. One possibility that has attracted more attention in recent years is that DM interacts with the SM via a massive spin-2 mediator [6][7][8][9][10][11][12][13][14][15] 1 . To date, there is no renormalizable model for the origin and nature of a massive spin-2 particle. Therefore, most studies either directly follow a simplified effective field theory approach or consider an origin from gravity in compactified extra-dimension, see however [11,12].
In these models, an important complication arises once external spin-2 particles are considered. Similar to the longitudinal modes of massive vector bosons the longitudinal part of the polarization vectors of massive spin-2 fields, µν 0 , exhibit an enhancement for large energies with µν 0 ∝ E 2 /m 2 where E is the energy and m the mass of the spin-2 particle G. In a theory with a single massive tensor field, the amplitude for the pair production of G from matter fields grows as E 6 in the high energy limit. This leads to a strong growth of the cross section ∝ E 10 which has a major impact on DM production in the early Universe if processes with √ s ≥ m G are relevant. Depending on the production mechanism the parameter range in which this is the case changes. On the one hand, for DM produced via thermal freeze-out, this high energy enhancement is important if its mass m DM exceeds the mass of the mediator, see e.g. [7,9]. On the other hand, for a freeze-in scenario, the relevant quantity is the temperature T instead and a strong increase of the production rate for T ≥ m G was found in e.g. [13,15].
Such a large growth of the amplitudes heralds a breakdown of the effective theory at energies well below the explicit scale Λ that enters in the effective operator. This is not necessarily problematic in a general effective theory and merely indicates that the unspecified UV completion has to become relevant at lower scales than naively expected. However, it is surprising in extra-dimensional models since they ought to remain perturbative up to the scale of higher dimensional gravity. These clashing expectations can be resolved by realizing that the 4D theory that arises from the KK-decomposition of extra dimensions does not just contain a single spin-2 field but a tower of KK-gravitons and their interactions instead. We recently analyzed the high energy limit of KK-graviton production in warped extradimensions and found that the contribution of the other KK-gravitons conspire to cancel the contributions that grow faster than E 2 [18] 2 . Naturally, this will have an impact on the DM production rates and can change the amount produced in the early Universe considerably. Therefore, it is of great interest to revisit the production of DM via spin-2 mediators in the context of extra-dimensions. For concreteness, we work with the warped extra-dimensions model of Randall and Sundrum [24]. We expect that our qualitative conclusions regarding the high energy behavior of the DM annihilation (or production) rates will also hold in extradimensional models with a different geometry. The results of our phenomenological study are clearly limited to the model we consider.
The structure of this paper is as follows. First, we introduce the Randall-Sundrum (RS) model which acts as our blueprint for spin-2 mediators from higher dimensional gravity. In Sec. 3 we discuss the production of DM in the early Universe focusing on the freeze-out mechanism. We pay particular attention to the modifications of the annihilation cross sections that arise once the full KK-tower is considered. Next, we reassess the phenomenology of thermally produced DM in the Randall-Sundrum model. Finally, we present our conclusions in Sec. 5. Additional material concerning the decay width of the graviton and radion and the derivation of some relations between the couplings of the KK-gravitons is presented in two Appendices.

The Randall-Sundrum model
The RS model consists of a 5D space compactified under a S 1 /Z 2 orbifold symmetry [24]. After compactification, one finds a 5D bulk space bounded by two branes. Only gravity propagates in the 5th dimension while SM fields are localized on one of the branes. To allow for a DM candidate, we extend the model minimally and add a single, only gravitationally interacting scalar that is localized on the same brane as the SM. In the following, we will briefly summarize the relevant ingredients of the RS model. A more pedagogical and detailed introduction can for example be found in [25][26][27].
The total action of the model is and includes a contribution from the bulk and the UV-and IR-branes, respectively. Only gravity propagates in the fifth dimension and, therefore, the bulk contribution reads 2 A related issue has also received attention in the context of spin-2 scattering. In a theory with a single massive graviton, the amplitude of GG → GG grows like s 5 . It has been shown that the tree-level exchange of any number of spin-0 and spin-1 fields cannot remove this behavior [19] but in extra-dimensional theories the growth is reduced to E 2 once the full KK-tower of massive gravitons is included, see [20] or [21][22][23].
where M 5 is the Planck mass in 5D while G is the determinant of the higher dimensional metric G M N and R denotes the Ricci scalar. We use large Latin letters, e.g. M, N , to indicate indices that run over all five dimensions while small Greek letters, e.g. µ, ν, only run over the usual four dimensions. The coordinate of the full 5D space, therefore, reads x M = (x µ , y). In the equation above, y has been replaced by a dimensionless coordinate ϕ = y/r c where r c is a measure of the size of the extra-dimension. To get the desired phenomenology, a bulk vacuum energy term Λ B has to be included. The contributions of the two branes located at y = 0 and y = πr c , respectively, are given by where g UV/IR are the metrics induced by G on the branes, L UV/IR indicate the corresponding matter fields Lagrangians and V UV/IR parametrize the branes' vacuum energies. We consider the case where both the SM and the DM are localized on the IR brane, i.e. L IR = L SM +L DM , while L UV = 0 for simplicity. To keep our model minimal we introduce a single scalar φ that only interacts gravitationally such that Solving Einstein's equation in the absence of matter and choosing the vacuum energy contributions such that 4D Poincaré invariance is realized on the branes leads to an infinitesimal line element where η µν is the Minkowski metric in 4D and the warping parameter k ≡ −Λ B 6 . It is often convenient to work with dimensionless quantities instead and one defines µ = kr c .
We want to work in an effective field theory in 4D which captures the low energy behavior of the model outlined above. Three steps are required to derive this EFT. First, one performs a weak-field expansion of the 5D metric in terms of an expansion parameter κ = 2/M 3/2 5 , that leads to a Lagrangian involving a tensor perturbationĥ µν , i.e. the 5D graviton, and a scalar perturbationr, the so-called radion. A detailed discussion of this expansion and the resulting Lagrangian up to third order in the fields can for example be found in [18,21]. Second, one utilizes a Kaluza-Klein decomposition of the 5D graviton to express it in terms of a tower of 4D spin-2 fields h (n) µν (x), which do not depend on the 5th dimension, and the wave functions ψ n (y) along the 5th dimension. The wave functions are fixed by a differential equation [28] 1 where m n is the mass of the n-th graviton and A(ϕ) = e −µ|ϕ| . A successful description of 4D gravity demands that the lightest graviton is massless and fixes the normalization of the zero mode ψ 0 = µ 1−e −2µπ . The other masses and the associated wave functions are obtained by solving eq. 2.7 with boundary conditions ∂ ϕ ψ n | ϕ=0,π = 0. The general solutions can be expressed in terms of Bessel-J and Bessel-Y functions and read (2.8) The normalization constants N n are fixed by the scalar product For e −µπ 1 the second term can be neglected and the solutions simplify to ψ n (ϕ) e 2µ|ϕ| N n J 2 γ n e µ(|ϕ|−π) and m n kγ n e −µπ , (2.10) where γ n denotes the nth zero of J 1 (x). For convenience it is conventional to perform a similar change in the normalization of the radionr (x) = 1 √ r c ψ r r(x) (2.11) even though ther field does not possess a y-dependence. Now we have a Lagrangian in which the fields h µν , r are already 4-dimensional but we still need to perform the integration over the 5th dimensions to arrive at a fully 4D theory. This last step fixes the coefficients of the interactions in terms of integrals over the 5D wave functions. The leading term in the interaction between gravitons and matter reads where T µν is the energy-momentum tensor of the matter field. Matching the massless graviton contribution to the equivalent term in GR fixes the relation between the (reduced) Planck mass in 4D, M P l and the parameters of the 5D theory Plugging in explicit expressions this leads to M 2 P l = M 3 5 k 1 − e −2µπ M 3 5 /k. The massive gravitons have a different normalization compared to the zero mode. This generates a different effective scale for these interaction which is given by Λ = M 3/2 5 √ r c /ψ n (π). In the large µ limit this reduces to Λ = M P l e −µπ such that the interaction strength with the massive gravitons is exponentially larger than with the massless one. In this limit, the radion contribution to the interaction Lagrangian is given by where T = η µν T µν is the trace of the energy-momentum tensor of the matter field. In addition, we need the three-point interactions between gravitons and radions and the next higher term in the expansion of the interaction with matter. The Lagrangian and the Feynman rules for these can be found in [18]. Here we only briefly comment on the coefficients of these interactions. The couplings between KK-gravitons and the radion can be obtained by plugging into the Lagrangian the solutions ψ n and ψ r and integrating out the 5D. In the case of interest only the 3-particles couplings are relevant, and in the limit e −µπ 1 they reduce to Λ −1 times a numerical coefficient that depends on the type of particles involved. They are summarized in Fig. 1. The numerical coefficients are given by The numerical coupling of the three radions vertex is trivially found to be χ rrr = 1/2.
In the original RS model, the radion is massless which implies a long-range force much stronger than gravity. This is clearly in contradiction with observations, and, therefore, a realistic model has to include a radion mass. One popular possibility to generate as mass is the Goldberger-Wise mechanism [29] which relies on a new bulk scalar to stabilizes the distance between the branes. Alternative ideas for the stabilization of extra-dimensions based on other kinds of bulk matter have also been put forwards, see e.g. [30,31]. We remain agnostic as to the origin of the radion mass and introduce its mass by hand.
Finally, it should be mentioned that a singlet scalar could also interact with the SM via the Higgs portal operator This interaction is very well studied, see e.g. [32,33], and we do not include it in our analysis since we want to keep our focus on the implications of higher dimensional gravity.

Dark Matter Production
Astrophysical and cosmological observations have clearly established the need for a large DM component in our Universe. The best measurement of its abundance comes from observations of the Cosmic Microwave Background (CMB) with the final analysis of the full Planck data yielding Ωh 2 = 0.1200(12) [1]. This result continues the best-known aspect of DM and, consequently, correctly reproducing it is one of the key marks of a successful theory. There are a number of mechanisms that explain how DM can be created in the early Universe. In the following, we consider the thermal freeze-out scenario while other possibilities such as a freeze-in will be discussed elsewhere. In this picture, DM is in thermal equilibrium with the SM plasma in the early Universe and decouples once the Hubble rate H becomes comparable to the creation and destruction rates of the DM particles in the thermal environment. The evolution of the DM number density n is described by the Boltzmann equation [34] n + 3Hn = − σv n 2 − n 2 eq , where σv is the thermally averaged cross section while n eq is the DM equilibrium density.
For temperature T small compared to the DM mass m, σv is given by [35] σv = 1 where K n (z) is the n-th Bessel-K function. In the second line, we have rewritten the expression in terms of the dimensionless variables w = s/4m 2 and x = m/T for convenience. Away from special kinematic configuration such as thresholds or resonances σv can often be Here X andX should be read as a placeholders for all kinds of SM particles.
approximated by a polynomial in v 2 or, equivalently, 1/x. Due to the KK-tower resonances and thresholds are ubiquitous in our model and we always use the full thermal average in our analysis. We also do not report approximate expressions in the velocity expansion but prefer to present σ(s) instead. Two classes of final states contribute to the φ annihilation rates, SM particles, and gravitational fields (i.e. gravitons and radions).

SM final states
Annihilations into SM fields are mediated by the exchange of KK-gravitons and the radion in the s-channel 3 . A set of representative Feynman diagrams is shown in Fig.2. If the DM is heavy such that the masses of SM particles can be neglected one finds a relatively simple expression for the cross section where Γ n is the width of the graviton. In this limit a good approximation of the width is Γ n ≈ 0.0968 m 3 n Λ 2 ; more detailed expressions that take all masses into account can be found in Appx. A. The whole KK-tower contributes and the cross section is characterized by a series of resonances. The thermal average smooths these out a bit but the resonance structure remains clearly visible in σv . In addition, there is a contribution from radion exchange which does not interfere with the gravitons. If the masses of the final state particles are neglected one finds where we have neglected the contribution of massless gauge bosons. These do not couple to the radion at tree level and their contribution at freeze-out is small. Expressions for the cross 3 In principle also the massless graviton contributes but since its interaction strength is fixed by M P l it can be neglected here. We fixed the temperature to T = m φ /20 which is in the ballpark of the typical freeze-out temperature.
section into different final states from both graviton and radion exchange with the full mass dependence can be found in [9]. For illustration, we show σv from radion and KK-graviton exchange for a representative set of parameters and two radion masses in Fig. 3. The radion cross section away from the resonance is dominantly s-wave, and, therefore depends only mildly on temperature. In contrast, the graviton cross section is either dominated by the resonances or has the leading contribution at d-wave, i.e. σv ∝ T 2 /m 2 φ . Therefore, the temperature has a very strong impact on the thermal average in this case. For illustration, we use T = m φ /20 which is in the ballpark of the expected freeze-out temperature. As can be seen, the radion and graviton cross section are comparable and the relative importance is set by the position of the resonances peaks. The radion peak is more pronounced than the graviton peaks since the radion width is smaller but also the impact of the different gravitons can be seen in the form of distinct peaks. At √ s much bigger than m 1 these start to wash out slightly and the cross section approaches a continuum that is dominated by graviton exchange.

Graviton and radion final states
The situation is more complicated when gravitons and radion in the final states are considered. In contrast to annihilation into SM fields which only depend on the interaction between the gravitons and the energy-momentum tensor we now need to include also the self-interaction of the gravitons, the radion-graviton coupling, and the effective four-point interactions φφG i G j , φφG i r and φφrr in our analysis. A set of representative Feynman diagrams for the case of annihilations into two gravitons is shown in Fig.4. The relevant Feynman rules and expressions for the coefficients of the graviton/radion interactions are reported in the appendix of [18]. Treating processes with external gravitons correctly is subtle. In the following, we briefly recapitulate the results of [18] where we discussed the challenges associated with graviton production in the Randall-Sundrum in detail, see also [20,21] for a related discussion on graviton scattering.
The longitudinal mode of the polarization tensor of massive spin-2 fields is proportional to s/m 2 which leads to a strong growth of the matrix elements with external gravitons in the high energy limit. Taking only the lightest graviton into account the matrix element φφ → G 1 G 1 grows as M ∝ s 3 . Specializing to the case of DM annihilations this would imply a very strong growth of the associated cross section with m φ , see for example [7,10]. This behavior is not altogether unexpected since an even worse growth of the amplitudes has been observed in spin-2 particle scattering amplitudes before [19,36,37]. In the absence of additional new physics, this leads to a breakdown of perturbative unitarity at a scale well below the naive effective scale of the theory Λ [38]. This is puzzling as the 5D graviton whose 4D decomposition leads to the Randall-Sundrum model does not have this issue. The seeming contradiction here can be resolved by realizing that the 4D theories derived from higher dimensional models consist of a KK-tower of gravitons and not an isolated massive spin-2 field. Once, the full content of the theory is considered the contributions of the different gravitons to s-channel exchange conspire to cancel the contributions to the matrix elements that grow faster than s [18,20,21]. This is conceptually similar to the unitarization of massive vector boson scattering in the SM [39] and has been dubbed unitarization by geometry in [20]. The chief difference with the usual Higgs mechanism is that a single field is not sufficient to restore perturbativity of the amplitudes in the EFT and the whole KK-tower is needed to fill the role of the SM Higgs 4 .
This has important implications for the annihilation rates into pairs of gravitons and graviton-radion final states. Unless the full tower of KK-gravitons is considered in the matrix elements unphysical contributions with a strong dependence on the center of mass energy arise. These can easily dominate the full cross section and lead to erroneous results. In practice, it is difficult to perform the sum over all internal KK-gravitons analytically. Therefore we use numerical methods when we want to include the full mass dependence. In the following, numerical results are calculated with a truncated KK-tower that includes the first 20 gravitons. This suppresses the prefactors of the unphysical contributions substantially and ensures that their contribution to the annihilation rates can be neglected throughout the parameter space considered here. Interestingly, in the high energy limit, analytic relations between the KK-graviton couplings derived in [18] can be applied at the cross section level. We find a simple expression for the cross-section for G n -pair production in this limit Thus, the high energy limit of σ(φφ −→ G n G n ) is universal and does not depend on the produced graviton. Analogous calculations for the production of different gravitons lead to (3.7) Similar considerations hold for the production of a radion and a KK-graviton. In the high energy limit the cross-section is found to be

8)
4 Additional complications arise if the scalar mode of the 5D graviton acquires a mass but it has recently been shown that these results hold if the popular Goldberg-Wise mechanism for generating the radion mass in the Randall-Sundrum model is treated carefully [40].
where the equality holds since the radion-graviton-graviton coefficients involved fulfill the relation (see Appx. B) Lastly, the cross-section for radion pair production is where the equality holds due to the fact that (see Appx. B) In the regime that is most relevant for freeze-out, i.e. √ s ≈ m φ , we cannot use the above approximations directly. However, one can instead parametrize s = 4m 2 φ w and perform an expansion in m φ . Keeping only the leading contributions in the small parameter w − 1 one arrives at Just as before the cross-sections are not sensitive to the type of KK-graviton(s) produced provided m φ is significantly bigger them m n . Furthermore, they are independent of the mass of the gravitons, i.e. of m 1 . These simple expressions can be easily thermally averaged using Eq. 3.2. For illustration we show σv for an exemplary set of parameters in Fig. 5. The full numerical result is depicted by the solid red line. To put this into perspective we also include a naive approximation that truncates the KK-tower after the heaviest graviton in the final state (black dashed). As can be seen, the naive results grow very strongly with m φ as expected from the high energy limit of the matrix elements discussed above 5 . In contrast, the full result exhibits a much more modest m φ dependence. For m φ ≈ 2m 1 the results already differ by an order of magnitude and, consequently, computations that ignore the sum over the KK-tower can not be trusted almost immediately after the channel becomes kinematically accessible. In addition, we have also included the analytic large m φ , small (w−1) approximation introduced In blue (dot-dashed) we show our approximation for the high m φ limit. The temperature is taken to be T = m φ /20 as in Fig. 3.
above (blue, dashed). While this result is not able the capture to behavior for masses close to the threshold, the agreement with the numerical computation is quite good for m φ 3m 1 and it reproduces the large mass behavior correctly.

Total cross section
Let us now turn to the total annihilation cross section. The cross section into SM particles is well described by Eq. 3.3 and Eq. 3.4 provided m φ m t . In the case of graviton and radion final states, the situation is a little bit more involved. Here, we have to take into account that more channels open as m φ increases. The analytic approximations do not capture the channels with m n close to m φ . Nevertheless, if more channels are open the bulk of them can be captured by Eqs. 3.12-3.15. and the total cross section in RS-particles is given by where G N is taken to be the heaviest graviton pair that can be produced. Due to their larger multiplicity, the sums are dominated by the mixed G i G j final. Identifying N = √ s/2∆m where ∆m is the mass spacing between the gravitons we find an estimate of total cross section into gravitational fields that captures the scaling in the large mass limit. The parameter c is an O(1) number that corrects for the fact that not all final state gravitons masses can be neglected in this sum. Weighting every channel by the phase space factor λ(s, m 2 i , m 2 j )/s, where λ is the Källén function, leads to the estimate c ≈ 0.6 for large N . In the more relevant regime of N 20 c is slightly smaller and increases with N . This approximation agrees well with our numerical results in the high mass regime.
In Fig. 6 we show a comparison of the thermally averaged annihilation cross section into SM and RS final states. As can be seen the cross sections into SM particles dominates the total cross section for the temperature shown here. At low masses, the RS cross section is dominated by radion pair production before the first graviton radion and the graviton pair production take over. In this mass range, the final state with the largest accessible mass dominates but already at the third peak several processes have to be added and it is no longer possible to identify a single leading channel. Throughout the whole mass range the SM cross section exceeds the one into gravitons and radions. In the high mass regime, where the difference is smallest, σv SM ≈ 3.5 × σv RS at x = 20 which decreases to σv SM ≈ 1.5 × σv RS at x = 30. As these ratios are taken in the high-energy regime, they are independent of m 1 and Λ. Note, however, that this is close to the breakdown of perturbative unitarity. Typically, freeze-out occurs at x ≈ 25. It is thus reasonable to approximate the full cross section by the SM contribution as long as m 1 m t and a O(30%) uncertainty in Ωh 2 is acceptable.

Phenomenology
In this section, we first review some of the most stringent experimental and theoretical limits on the parameters of our model. Then we combine these with results for the cosmologically preferred parameter space of thermally produced DM and identify the regions that allow for a successful thermal freeze-out.

Collider limits
A stringent bound on the RS model comes from LHC experiments, in particular from searches for a heavy di-photon resonance. We use the limits based on 37fb −1 of √ s = 13 TeV data reported by ATLAS [41]. In principle one could also consider the limits on the radion from searches for a spin-0 resonance. In practice, we find that these are not as strong and omit them from our analysis.
To extract limits on the parameters of the theory we compute the RS di-photon cross section with CalcHEP [42] and compare it to the limit presented by the collaboration. The experimental results are reported for some fixed value of k/M P l which impacts the width of the resonance. However, throughout most of the parameter space the width of the graviton is small compared to the energy resolution and in the region where this is no longer the case the limit is insensitive to the width. Therefore, we can use the limits of [41] directly. . The dashed lines correspond to the production of radion-radion (purple), G 1 -radion (blue), G 1 G 1 (green) and G 1 G 2 (orange). The contribution from heavier gravitons is not shown separately but included in the black line. The temperature is taken to be T = m φ /20 as in Fig. 3. For illustration we show the limit on Λ as a function of m 1 in Fig. 7. As can be seen, the bound is particularly stringent at low graviton masses and exceeds 100 TeV for m 1 < 1 TeV. It stays above 50 TeV up to m 1 ≈ 2 TeV. Once the mass of the graviton exits the energy window accessible for a resonance search, i.e. for m 1 ≥ 2.5 TeV, the limit softens considerable. For m 1 ≥ 5 TeV the collaboration does not report any limits in the Randall-Sundrum model. In principle, one could attempt to extend them by recasting a limit on virtual graviton exchange in large extra-dimensions that comes from a different analysis of a subset of the full data. We do not attempt this here since we will find that the part of the parameter space where this becomes relevant is disfavored by theoretical considerations.

Perturbative unitarity
Since we are working in an effective theory it is of great importance to understand where our calculations are valid. The proportionality of the vertices' couplings between matter and gravitons (and between gravitons and radions) to Λ −1 suggests that perturbativity ought to break at this scale. However, the KK decomposition obscures the fact that the 5D theory . from which our model arises is already an effective theory with a breakdown scale that is set by 5D gravity with M 5 and the warp parameter k as dimensionful parameters. Therefore, one might wonder which scale actually controls the validity of the 4D theory. One possibility for assessing this question is perturbative unitarity. We use partial wave perturbativity following the approach of e.g. [39,43,44]. The helicity matrix element for the J-th partial wave with initial and final states i and f is given by where β if denotes a kinematical factor, θ is the scattering angle while d J µµ (θ) is the J-th Wigner-functions for total spins of initial and final states µ and µ . Unitarity of the S-matrix enforces an inequality on the elastic elements which can be turned in an equality by diagonalizing the matrix M J ij . Therefore, the limit inferred from a single elastic scattering channel is conservative and a more stringent bound can be obtained by considering the full set of accessible states. For simplicity we focus on the J = 0 amplitude with µ = µ = 0 and start with the elastic scattering of the scalar field φ, see Fig. 8 for a representative set of the Feynman diagrams. For the sake of simplicity, we focus on the s-channel diagram for now 6 . The contribution from the full KK-tower to this amplitude is given by Naively, one could try to estimate the sum by assuming that the first N -gravitons have negligible mass compared to √ s, i.e. m n √ s. Taking Γ n m n and neglecting the contribution from heavier gravitons suggest the following scaling The number N is approximately given by the total energy divided by the mass-gap between each graviton ∆m, which tends to a constant as the number of the related gravitons grow, i.e.
It follows that the sum of the propagators does not contribute as 1/s in the high energy regime, but rather as N/s ≈ 1/ √ s, which makes the amplitude grow ∝ s 3/2 . A more rigorous way to arrive at this conclusion is to consider the large s approximation of S(s) [45,46] which confirms our rough estimate. In contrast to the hand-waving argument above, the validity of this approximation is easily quantified and found to agree to better than 5% for This implies a scaling M ∝ s 3/2 which points towards a breakdown of unitarity at a scale below Λ. Using the approximation eq. 4.6 and including the t-and u-channel contribution Figure 9: Largest eigenvalue of the graviton scattering matrix λ max rescaled by Λ 2 /s as a function of N . The red-dashed line shows the linear fit that was used to extract the perturbativity limit.
yields a perturbativity limit of Thus, the breakdown scale is set by the warped 5D Planck massM 5 . A similar argument also holds for DM annihilations into SM particles. In this case the growth of M is also reflected in the large √ s behavior of the σ SM ∝ s 2 , or, equivalently σv SM ∝ m 4 φ , that can be see in the upper right corner of Fig. 6.
Even though it is intuitive that the scale up to which the 4D theory remains valid is set by the scale in the underlying 5D theory, it is slightly disconcerting that the result relies on the high energy limit of the sum of graviton propagators that only becomes a good approximation at √ s barely in the perturbative regime. Therefore, we aim to strengthen this result by considering a different part of the S-matrix, namely scattering of KK-gravitons. In the case of G i G i → G i G i most of the work has already been done by [21]. As mentioned above the leading contribution to this matrix element is proportional to s 1 which might cast some mild doubt on the s 3/2 scaling found in φ scattering. However, one has to take into account that increasing √ s also opens the thresholds for the production of more graviton pairs. Thus we cannot neglect inelastic scattering and we perform a coupled channel analysis following the approach of [43]. We combine the results given in [21] with our own computation to derive the inelastic scattering matrix elements for G i G j → G k G l . As noted previously, these matrix elements grow faster than s unless the appropriated sum rules for the coefficients of the three-and four-graviton are applied. Unfortunately, we have not found a simple analytic expression for these so we opt for a numerical analysis instead. We confirm that the scaling of the amplitudes is reduced to O(s) if a sufficient number of exchanged-gravitons are included. For computational reasons and since it suffices to our purpose, we consider only the 0-helicities of the gravitons and construct a matrix with all states whose production is allowed at √ s. Numerical evaluation of the partial wave amplitudes then leads to the following scaling for the largest eigenvalue for √ s m 1 , where N is the number of the heaviest graviton that can be pair produced. An illustration of this behavior together with our numerical results as can be seen in Fig. 9. Recalling that N ≈ √ s/2∆m and extracting the proportionality constant from the numerical results we find the perturbativity condition to be √ s 1.9M 5 . (4.10) This confirms and improves the previous result and, more importantly, is independent of the high energy limit of the sum of graviton propagators.

Parameter space favored by freeze-out
We focus our analysis on DM and lightest KK-graviton heavier than SM particles since this region potentially allows for the production of DM via the freeze-out mechanism. It would be interesting to investigate to which extent a lighter and more weakly coupled lightest graviton is in agreement with experiments but since this is not expected to be consistent with thermally produced DM we do not consider this case further. The annihilation cross section in this regime is dominated by SM final states, see Sec. 3.3. For the efficient computation of the relic density, we implement the annihilations into SM particles via the first six gravitons in MicrOMEGAs [47,48] and remove the artificially large contribution from RS final states by hand 7 . We scan over m φ and m 1 while keeping the radion mass fixed to 1 TeV and extract the value of Λ required for the thermal production of the observed relic density. The results are visualized in Fig. 10. The dashed contours indicate lines of constant Λ. On top of this, we superimpose the limits from collider searches and theoretical consistency conditions. In the blue region starting at m φ 2 TeV, the value of Λ preferred by freeze-out is in the non-perturbative region and hence our computation is not reliable. On the other hand, the strong collider bounds of Fig. 7 excludes the thermal DM region with m 1 ≤ 5 TeV almost completely. The only exception is for m φ ≈ m r /2 where the radion resonance allows for efficient annihilations even for comparatively large values of Λ. In the top half, the majority of the parameter space is characterized by Γ 1 ≥ m 1 /2, which also casts severe doubts on the The blue (orange) region is outside of the range of validity of perturbative computations due to breakdown of partial wave unitarity (large graviton width, γ 1 ≥ m 1 /2). In the white region thermal production of the observed relic density is viable. reliability of a perturbative computation. This leaves two regions where thermal production remains viable, i.e. close to the radion resonance, and for m 1 slightly above 5 TeV. This second small corner is uncomfortably close to the theoretical consistency conditions. While we can not exclude that a more complicated underlying theoretical model justifies 2m φ = m r or that the perturbativity limits are close to being violated in nature, this makes it unlikely that thermal production in the minimal Randal-Sundrum is the origin of DM. Naturally, this statement cannot be generalized to extra-dimensional theories as a whole and more complex extra-dimensional constructions with fields beyond the graviton in the bulk remain viable, see for example [8,49].

Conclusions
There is a great interest in new ideas for the interactions between DM and the SM. In this paper, we studied DM interacting with the SM via a massive spin-2 mediator. In this context, an interesting complication arises as soon as the energy is sufficient to produce the mediators on-shell. The longitudinal modes of the polarization tensor of massive spin-2 fields are proportional to E 2 /m 2 which leads to the rapid growth of their production rate with energy almost instantly above the threshold. A similar behavior is known from studies of spin-2 scattering in massive gravity and is a sign of an early breakdown of perturbative unitarity in theories with a single massive tensor field. Due to the well-known issues of theories with spin-2 fields, all studies of these interactions have to resort to an EFT approach. Here, we focused on an EFT that is generated by the compactification of a warped extra-dimensional model. Interestingly, such a model is expected to remain perturbative up to the scale of 5D gravity. Relying on our earlier work on the unitarization of spin-2 production amplitude in warped extra-dimensions [18], we are able to show that the strong growth of the cross section above the threshold is not physical and gets reduced considerably once the full KK-tower is included in the computation. Thus, an EFT with a limited number of degrees of freedom is not able to capture the behavior of extra-dimensional theories, and the results obtained with this approximation are misleading. We derived the DM annihilation cross sections carefully taking the KK-tower into account and evaluated it numerically. To better analyze the high mass behavior and to cross-check our numerical study we also report approximate analytic expressions in this limit.
We have used our results for the annihilations to reassess the parameter space that allows for a thermal production of the observed relic abundance. After confronting these results with limits from collider searches for the lightest KK-graviton and theoretical consistency conditions, we find that gravitationally interacting scalar DM is under more pressure than previously anticipated. The only region which is clearly viable features resonant annihilations via the radion and is therefore characterized by an unexpected tuning between the radion and the DM mass. In addition, we find a second small region that is dangerously close to both the experimental limit and the theoretical consistency conditions. Given that a perturbative computation will likely be modified even before the breakdown of unitarity, it is difficult to make a definitive statement about this region without specifying the UV completion of extradimensional gravity. In any case, the upcoming LHC runs will improve the experimental limits and close the gap.
In this work, we mainly focused on the non-relativistic limit of the annihilation cross section which is relevant for thermal freeze-out. However, it is clear that the unitarization of the annihilation cross section by the KK-tower will also affect the cross sections in the regime that matters for freeze-in. A study of the impact of this production mechanism would therefore be very interesting. Finally, we also want to note that there is an increasing interest in theories with DM of spin-2 or higher [50,51]. While our analysis does not directly address these ideas one might expect that variations of the unitarity issues found here also persist there. Thus, the results derived with this approach could be more sensitive to unspecified UV physics than naively expected.

Acknowledgments
We thank the Max-Planck Institute of Physics (MPP) for access to some of the computational resources employed in the work.

Appendices A Decay widths
In this appendix, we present the decay widths of the gravitons and the radion. In both cases the decay into SM-particles and into other RS-particles in the final state are considered.

A.1 Graviton width
The graviton decays into all SM particles, into lighter gravitons, and radions. Our results for the partial width into SM particles match those reported in [10]. The width of a heavy KKgraviton decaying to a light KK-graviton in a different geometry has been reported in [52]. Even though their result ought to match ours once the coefficients of the KK-graviton vertex are adapted to our geometry, we do not find agreement. The difference has been traced to the relative sign of the part R 4 (g) and the contribution with a 5D derivative to the three-graviton vertex. We have checked our result by comparing it with the three-graviton vertex reported by [20]. In addition, we verified that the unitarization of the graviton production amplitude does not go through with the other sign.
The partial width for decays into a pair of Higgs bosons reads where N c is the number of colors. For the massive gauge bosons, Z andW ± the width is while the expression for the massless gauge bosons γ and g reads If the masses of the SM particles can be neglected this leads to combined width of The general expression for the decay into two lighter gravitons G m and G k is As long as n is not very large these contributions are small compare to the width into SM particle and can be neglected.

A.2 Radion Decay Widths
The partial width into Higgs boson is given by In the case of the massive gauge bosons Z and W ± the width is while the decay into massless gauge bosons γ and g is given by The radion does not couple to massless particles at tree-level. The loop contribution is encoded in the coefficients C EM and C S . These channels are not relevant to our analysis. A discussion of the coefficients can be found in [53]. In principle, the radion can also decay to gravitons if it is heavy enough. In this work, the radion is assumed to be lighter than the first graviton and we do not include these decay channels. to simplify the expressions of the cross section. In the following we provide an explicit proof for the interested reader.

B Coupling relations
As detailed in [18], c n,α and d n,α can be evaluated with the help of the Fourier-Bessel expansion, see for example [54]. If a function f (x) is continuous on [0, 1] such that f (1) = 0, the integral where γ ν,k is the k−th root of J ν and a ν,k = 2 J ν+1 (γ ν,k ) 2 1 0 du uf (u)J ν (γ ν,k u) . (B.7) Since we will work only with the roots of J 1 , we define γ k as the k−th root of J 1 and we will set also ν = 1 such that Our strategy consists of picking an appropriate function f (x), such that we can obtain the relations of interest. We will also use some relations derived in [18] and the definitions of the couplings given in eq. 2.15. In [55], it is proved that where γ ν,n is the nth-zero of the Bessel-J ν function. In the case of interest, ν = 1, this directly leads to e 0 = ∞ n=1 χ nrr = 1.