Gravity-mediated scalar Dark Matter in warped extra-dimensions

We revisit the case of scalar Dark Matter interacting just gravitationally with the Standard Model (SM) particles in an extra-dimensional Randall-Sundrum scenario. We assume that both, the Dark Matter and the Standard Model, are localized in the TeV brane and only interact via gravitational mediators, namely the graviton Kaluza-Klein modes and the radion. We analyze in detail the dark matter annihilation channel into two on-shell KK-gravitons, and contrary to previous studies which overlooked this process, we find that it is possible to obtain the correct relic abundance for dark matter masses in the range [1, 10] TeV even after taking into account the strong bounds from LHC Run II. We also consider the impact of the radion contribution (virtual exchange leading to SM final states as well as on-shell production), which does not significantly change our results. Quite interestingly, a sizeable part of the currently allowed parameter space could be tested by LHC Run III and by the High-Luminosity LHC.

The Nature of Dark Matter (DM) is one of the long-standing puzzles that still have to be explained in order to claim that we have a "complete" picture of the Universe. On the one side, both from astrophysical and cosmological data (see, e.g., ref. [1] and refs. therein), rather clear indications regarding the existence of some kind of matter that gravitates but that does not interact with other particles by any other detectable mean can be gathered.
On the other hand, no candidate to fill the rôle of Dark Matter has yet been observed in high-energy experiments at colliders, nor is present in the Standard Model (SM) spectrum. Within SM particles, the only ones that share with Dark Matter the property of being weakly coupled to SM matter are neutrinos. However, experimental searches have shown that neutrinos constitute just a tiny fraction of what is called non-baryonic matter in the Universe energy budget [2]. Most of the suggestions for physics Beyond the Standard Model (BSM), therefore, include one or several possible candidates to be the Dark Matter. Under the assumptions of the "WIMP paradigm" (with "WIMP" standing for "weakly interacting massive particle"), these new particles have in common to be rather heavy and with very weak interactions with SM particles. Two examples of these are the neutralino in supersymmetric extensions of the SM [3] or the lightest Kaluza-Klein particle in Universal Extra-Dimensions [4,5]. Searches for these heavy particles at the LHC have pushed bounds on the masses of the candidates to the TeV range, a region of the parameter space rather difficult to test for experiments searching for Dark Matter particles interacting directly within the detector (see, e.g., ref. [6]) or looking at annihilation products of Dark Matter particles [7]. Both for this reason and for the fact that very heavy WIMP's are relatively unnatural in theories that want to solve the hierarchy problem and not only host some Dark Matter candidates, models in which the Dark Matter particles are either "feebly interacting massive particles" (FIMP's) [8] or "axion-like" very light particles (see, e.g., ref. [9]) have been constructed. As a result, at present a very rich (and complicated) landscape of models explaining the Nature of Dark Matter exists, and experimental searches have to look for very different signals.
In this paper we want to explore in some detail a possibility that was advanced in the literature several times in the last ten to twenty years. The idea is that the interaction between Dark Matter particles and the SM ones, though only gravitational, may be enhanced due to the fact that gravity feels more than the standard 3 + 1 space-time dimensions. Extra-dimensional models have been proposed to solve the hierarchy problem, related to the large hierarchy existing between the electro-weak scale, Λ EW ∼ 250 GeV, and the Planck scale, M P ∼ 10 19 GeV. In all these models, the gravitational interaction strength is generically enhanced with respect to the standard picture since the "true" scale of gravitation is not given by M P but, rather, by some fundamental scale M D (where D is the number of dimensions). The two scales, M P and M D are connected by some relation that takes into account the geometry of space-time. In so-called Large Extra-Dimensions models (LED) [10][11][12][13][14], for example, warped extra-dimensions (also called Randall-Sundrum models) [15,16], on the other hand, the separation between M P and M D is not very large, M 2 P = 8π(M 3 D /k) [1 − exp(−2πr c k)], where k is the curvature of the space-time along the extra-dimension and r c is the distance between two points in the extra-dimension. However, all physical masses have an exponential suppression with respect to M P due to the curvature k, m = exp(−2πr c k) m 0 . In this picture, m 0 is a fundamental mass parameter of order M P and m is the mass tested by a 4-dimensional observer. In the ClockWork/Linear Dilaton model (CW/LD), eventually, the relation between M P and M D is a combination of a volume factor, as for LED models, and a curvature factor, as for warped models [17].
The possibility that Dark Matter particles, whatever they be, have an enhanced gravitational interaction with SM particles have been studied mainly in the context of warped extra-dimensions. The idea was first advanced in refs. [18,19] and subsequently studied in refs. [20][21][22][23][24][25]. As already stressed, the Nature of Dark Matter is still unknown. In particular, if new particles are added to the SM spectrum to act as Dark Matter particles, their spin is completely undetermined. In the publications above, therefore, scalar, fermion and vector DM particles have been usually considered. In this paper, on the other hand, we only consider scalar Dark Matter. We have been led to this decision by the fact that, maybe unexpectedly, we have found significant regions of the model parameter space for which the thermal relic abundance can be achieved and that can avoid present experimental bounds and theoretical constraints (in contrast, for example, with the conclusions of ref. [21]). Interestingly enough, most of the allowed parameter space will be tested by the Run III at the LHC and by its high-luminosity version, the HL-LHC. On the way to achieve the correct relic abundance, we have found some discrepancies with existing literature on the subject when looking for DM annihilation into Kaluza-Klein gravitons. In addition, in order to give a consistent picture of this possibility in the framework of warped extra-dimensions, we have also taken into account the DM annihilation through and to radions within the Goldberger-Wise approach [26], finding that this channel may also give the correct relic abundance, though in a very tiny region of the parameter space difficult to test at the LHC.
In forthcoming publications we plan to extend our study to DM particles with a different spin and explore other extra-dimensional scenarios, such as LED and CW/LD. The paper is organized as follows: in section 2 we outline the theoretical framework, reminding shortly the basic ingredients of warped extra-dimensional scenarios and of how dark matter can be included within this hypothesis; in section 3 we show our results for the annihilation cross-sections of scalar DM particles into SM particles, KK-gravitons and/or radions; in section 4 we review the present experimental bounds on the Kaluza-Klein graviton mass from LEP and LHC, as well as on the DM mass from direct and indirect search experiments, and we remind the theoretical constraints coming from unitarity violation and effective field theory consistency; in section 5 we explore the allowed parameter space such that the correct relic abundance is achieved for scalar DM particles; and, eventually, in section 6 we conclude. In the appendices we give some of the mathematical expressions used in the paper: appendix A contains the KK-graviton propagator and polarization tensor; in appendix B we provide the Feynman rules for our model; in appendix C we give the -2 -expressions for the decay amplitudes of the KK-graviton and of the radion; and, eventually, in appendix D we give the formulae for the annihilation cross-sections of dark matter particles into Standard Model particles, KK-gravitons and radions.

Theoretical framework
In this section, we shortly review the Warped Extra-Dimensions scenario (also called Randall-Sundrum model [15]) and introduce our setup to include Dark Matter in the model, we give the relevant formulae to compute the DM relic abundance and eventually provide the DM annihilation cross-sections into SM particles, Kaluza-Klein gravitons and into radions.

A short summary on warped extra-dimensions
The popular Randall-Sundrum scenario (from now on RS or RS1 [15], to be distinguished from the scenario called RS2 [16]) consider a non-factorizable 5-dimensional metric in the form: where σ = kr c |y| and the signature of the metric is (+, −, −, −, −). In this scenario, k is the curvature along the 5th-dimension and it is O (M P ). The length-scale r c , on the other hand, is related to the size of the extra-dimension: we only consider a slice of the space-time between two branes located conventionally at the two fixed-points of an orbifold, y = 0 (the so-called UV-brane) and y = π (the IR-brane). The 5-dimensional space-time is a slice of AdS 5 and the exponential factor that multiplies the M 4 Minkowski 4-dimensional spacetime is called the "warp factor". Notice that, in order to have gravity in 4-dimensions, in general η µν → g µν , with g µν the 4-dimensional induced metric on the brane. The action in 5D is: with M 5 the fundamental gravitational scale, G M N and R (5) the 5-dimensional metric and Ricci scalar, respectively, and Λ 5 the 5-dimensional cosmological constant. As usual, we consider capital latin indices M, N to run over the 5 dimensions and greek indices µ, ν only over 4 dimensions. The Planck mass is related to the fundamental scale M 5 as: whereM P = M P / √ 8π = 2.435 × 10 18 GeV is the reduced Planck mass. We consider for the two brane actions the following expressions: where f IR , f UV are the brane tensions for the two branes. In Randall-Sundrum scenarios, in order to achieve the metric in eq. (2.1) as a classical solution of the Einstein equations, the brane-tension terms in S UV and S IR are chosen such as to cancel the 5-dimensional cosmological constant, f 4 IR = −f 4 UV = −24M 3 5 Λ 5 . Throughout this paper, we consider all the SM and DM fields localized on the IR-brane, whereas on the UV-brane we could have any other physics that is Planck-suppressed. We assume that DM particles only interact with the SM particles gravitationally and, for simplicity, we focus on scalar DM. More complicated DM spectra (with particles of spin higher than zero or with several particles) will not be studied here. Notice that, in 4dimensions, the gravitational interactions would be enormously suppressed by powers of the Planck mass. However, in an extra-dimensional scenario, the gravitational interaction is actually enhanced: on the IR-brane, in fact, the effective gravitational coupling is Λ = M P exp (−kπr c ), due to the rescaling factor √ G (5) / −g (4) . It is easy to see that Λ M P even for moderate choices of σ. In particular, for σ = kr c ∼ 10 the RS scenario can address the hierarchy problem. In general, we will work with Λ = O(1 TeV) but not necessarily as low as to solve the hierarchy problem. Expanding the 4-dimensional component of the metric at first order about its static solution, we have: with κ 5 = 2M −2/3 5 . The 5-dimensional field h µν can be written as a Kaluza-Klein tower of 4-dimensional fields as follows: The h n µν (x) can be interpreted as the KK-excitations of the 4-dimensional graviton. The χ n (y) factors are the wavefunctions of the KK-gravitons along the extra-dimension. Notice that in the 4-dimensional decomposition of a 5-dimensional metric, two other fields are generally present: the graviphoton, G (5) µ5 , and the graviscalar G (5) 55 . It has been shown elsewhere [15] that graviphotons are massive due to the breaking of 5-dimensional translational invariance induced by the presence of the branes. On the other hand, the graviscalar field is relevant to stabilize the size of the extra-dimension and it will be discussed below when introducing the radion.
The equation of motion for the n-th KK-mode is given by: where m n is its mass. Using the Einstein equations we obtain [27]: being J 2 and Y 2 Bessel functions of order 2 and z n (y) = m n /ke σ(y) . The N n factor is the n-th KK-mode wavefunction normalization. In the limit m n /k 1 and e kπrc 1, the coefficient α n becomes α n ∼ x 2 n exp (−2kπr c ), where x n are the are the roots of the Bessel function J 1 , J 1 (x n ) = 0, and the masses of the KK-graviton modes are given by: (2.12) Notice that, for low n, the KK-graviton masses are not equally spaced, as they are proportional to the roots of the Bessel function J 1 . This is very different from both the LED and the CW/LD scenarios, however for large n the spacing between KK-graviton modes become so small that all extra-dimensional scenarios eventually coincide, m n ∼ n/R (being R some relevant length scale specific to each scenario). The normalization factors can be computed imposing that: In the same approximation as above, i.e. for m n /k 1 and e kπrc 1, we get: (2.14) Notice the difference between the n = 0 mode and the n > 0 modes: for n = 0, the wave-function at the IR-brane location y = π takes the form whereas for n > 0: The important difference can be easily understood by looking at the coupling between the energy-momentum tensor and gravity at the location of the IR-brane: where the only scale is the fundamental gravitational scale M 5 . However, if we separate the n = 0 and the n > 0 modes we get: from which is clear that the coupling between KK-graviton modes with n = 0 is suppressed by the effective scale Λ and not by the Planck scale.

JHEP01(2020)161
It is useful to remind here the explicit form of the energy-momentum tensor: where we have introduced the scalar singlet field S to represent the DM particle in our scenario.
Notice that a scalar DM particle will also interact with the SM through the so-called "Higgs portal", namely since this term is always allowed. However, such coupling is strongly constrained (see section 2.3), and we neglect its effect in our analysis.

Adding the radion
Stabilizing the size of the extra-dimension to be y = πr c is a complicated task. In general (see, e.g., refs. [28][29][30]) bosonic quantum loops have a net effect on the border of the extradimension such that the extra-dimension itself should shrink to a point. This feature, in a flat extra-dimension, can only be compensated by fermionic quantum loops and, usually, some supersymmetric framework is invoked to stabilize the radius of the extra-dimension (see, e.g., ref. [31]). In Randall-Sundrum scenarios, on the other hand, a new mechanism has been considered: if we add a bulk scalar field Φ with a scalar potential V (Φ) and some ad hoc localized potential terms, δ(y = 0)V UV (Φ) and δ(y = πr c )V IR (Φ), it is possible to generate an effective potential V (ϕ) for the four-dimensional field ϕ = f exp (−kπT ) (with f = 24M 3 5 /k and T = r c ). The minimum of this potential can yield the desired value of kr c without extreme fine-tuning of the parameters [26,32].
As in the spectrum of the theory there is already a scalar field, the graviscalar G 55 , the Φ field will generically mix with it. The KK-tower of the graviscalar is absent from the low-energy spectrum, as they are eaten by the KK-tower of graviphotons to get a mass (due to the spontaneous breaking of translational invariance caused by the presence of one or more branes). On the other hand, the KK-tower of the field Φ is present, but heavy (see ref. [33]). The only light field present in the spectrum is a combination of the graviscalar zero-mode and the Φ zero-mode. This field is usually called the radion, r. Its mass can be obtained from the effective potential V (ϕ) and is given by where v v is the value of Φ at the visible brane and = m 2 /4k 2 (with m the mass of the -6 -JHEP01(2020)161 field Φ). Quite generally 1 and, therefore, the mass of the radion can be much smaller than the first KK-graviton mass.
The radion, as for the KK-graviton case, interacts with both the DM and SM particles. It couples with matter through the trace of the energy-momentum tensor T [18]. Massless gauge fields do not contribute to the trace of the energy-momentum tensor, but effective couplings are generated from two different sources: quarks and W boson loops and the trace anomaly [34]. Thus the radion Lagrangian takes the following form [33,35]: where F µν , F a µν are the Maxwell and SU c (3) Yang-Mills tensors, respectively. The C 3 and C EM constants encode all information about the massless gauge boson contributions and are given in appendix B.

The DM relic abundance in the freeze-out scenario
Experimental data ranging from astrophysical to cosmological scales point out that a significant fraction of the Universe energy appears in the form of a non-baryonic (i.e. electromagnetically inert) matter. This component of the Universe energy density is called Dark Matter and, in the cosmological "standard model", the ΛCDM, it is usually assumed to be represented by stable (or long-lived) heavy particles (i.e. non-relativistic, or "cold"). Within the thermal freeze-out scenario the DM component is supposed to be in thermal equilibrium with the rest of particles in the Early Universe. The evolution of the dark matter number density n DM in this paradigm is governed by the Boltzmann equation [36]: where T is the temperature, H(T ) is the Hubble parameter as a function of the temperature, and n eq DM is the DM number density at equilibrium (see, e.g., ref. [36]). The Boltzmann equation is governed by two factors: one proportional to H(T ) and the second to the thermally-averaged cross-section, σv . In order for n DM (T ) to freeze-out, as the Universe expanded and cooled down the thermally-averaged annihilation cross-section σv times the number density should fall below H(T ). At that moment, DM decouples from the rest of particles leaving an approximately constant number density in the co-moving frame, called relic abundance.
The experimental value of the relic abundance can be computed starting from the DM density in the ΛCDM model. From ref. [37] we have Ω CDM h 2 = 0.1198 ± 0.0012, where h parametrizes the present Hubble parameter. Solving eq. (2.23), it can be found the thermally-averaged cross-section at freeze-out σ FO v = 2.2×10 −26 cm 3 /s [38]. Notice that for m DM > 10 GeV, the relic abundance is insensitive to the value of m DM and therefore the thermally-averaged annihilation cross section σ FO needed to obtain the correct relic abundance is not a function of the DM particle mass.
When comparing the prediction of a given model to the expectation in the freeze-out scenario, the key parameter to compute the relic abundance is, thus, σv . In order to -7 -JHEP01(2020)161 obtain this quantity, we must first calculate the total annihilation cross-section of the DM particles (represented in our case by the field S): (2.24) where in the first term, σ ve ("ve" stands for "virtual exchange"), we sum over all SM particles. The second term, σ GG , corresponds to DM annihilation into a pair of KK-gravitons, G n G m . Eventually, the third term, σ rr , corresponds to DM annihilation into radions.
If the DM mass m S is smaller than the mass of the first KK-graviton and of the radion, only the first channel exists. Since in the freeze-out paradigm the DM particles have small relative velocity v when the freeze-out occurs, it is useful to approximate the c.o.m. energy s as s ∼ 4m 2 S , and keep only the leading order in the so-called velocity expansion. Formulae for the DM annihilation into SM particles within this approximation are given in appendix D.
Notice that DM annihilation to SM particles can occur through three possible mediators: the Higgs boson, the KK-tower of gravitons and the radion. The first option, that depends on the coupling introduced in eq. (2.21), has been extensively studied. Current bounds (see for instance [39,40] for recent analyses) rule out DM masses m S 500 GeV (except for the Higgs-funnel region, m S m h /2) and future direct detection experiments such as LZ [41] will either find DM or exclude larger masses, up to O(TeV). In the presence of other annihilation channels, as in our case, if LZ does not get any positive signal of DM it will lead to a stringent limit on the Higgs portal coupling λ hS , so that the Higgs boson contribution to DM annihilation into SM particles will be negligible for DM masses at the TeV scale [40,42]. In the rest of the paper, we will assume that λ hS is small enough so as to be irrelevant in our analysis, and we will not consider this channel any further.
On the other hand, depending on the particular values for the radion mass (determined by the specific features of the bulk and localized scalar potentials) and the KK-graviton masses (fixed by k, M 5 and r c ), radion or KK-graviton exchange can dominate the annihilation amplitude. When computing the contribution of the radion and KK-graviton exchange to the DM annihilation cross-section into SM particles, it is of the uttermost importance to take into account properly the decay width of the radion and of the KKgravitons, respectively. 1 Notice that the DM annihilation cross-section into SM particles via virtual exchange of KK gravitons is velocity suppressed (d wave), due to the spin 2 of the mediators, while the corresponding one through virtual radion is s wave.
Within the Goldberger-Wise stabilization mechanism, the radion is expected to be lighter than the first KK-graviton mode, so the next channel to open is usually the DM annihilation into radions. The analytic expression for σ rr (S S → r r) in the approximation s ∼ 4m 2 S is given in appendix D. It is also s wave. Eventually, for DM masses larger than the mass of the first KK-graviton mode, annihilation of DM particles into KK-gravitons becomes possible and the last channel in JHEP01(2020)161 eq. (2.24) opens. As the KK-number is not conserved due to the presence of the branes in the extra-dimension (that breaks explicitly momentum conservation in the 5th-dimension), any combination of KK-graviton modes is possible when kinematically allowed. Therefore, we must sum over all the modes as long as the condition 2m S ≥ m Gn + m Gm is fulfilled. The analytic expression for σ GG (S S → G n G m ) at leading order in the velocity expansion is also given in appendix D, and it turns out to be s wave as well. Notice that we will not take into account annihilation into zero-modes gravitons, G 0 G 0 or G 0 G n , as these channels are Planck-suppressed with respect to the production of a pair of massive KK-graviton modes, G n G m .
As the velocity expansion approximation may fail in the neighbourghood of resonances and, in the RS model, the virtual graviton exchange cross-section is indeed the result of an infinite sum of KK-graviton modes, we computed the analytical value of σv using the exact expression from ref. [43]: where K 1 and K 2 are the modified Bessel functions and v is the relative (Møller) velocity of the DM particles.

Scalar DM annihilation cross-section in RS
For relatively low DM mass the only open annihilation channel is into SM particles through KK-graviton or radion exchange. Direct production of radions or KK-gravitons in the final state becomes allowed for DM mass m S ≥ m G 1 , m r , where m S and m G 1 are the DM and the first KK-graviton masses, respectively.

Virtual KK-graviton exchange and on-shell KK-graviton production
We plot in figure 1 the different KK-graviton contributions to σv separately, so as to understand clearly the main features. We consider first the case of DM annihilation into SM particles through KK-graviton exchange, summed over all virtual KK-gravitons, σ ve,G . This result is shown by the solid (purple) line in figure 1 as a function of the DM mass m S , for the particular choice Λ = 100 TeV and m G 1 = 1 TeV. When the DM particle mass is nearly half of one of the KKgraviton masses, s = 4m 2 S ∼ m 2 Gn , the resonant contribution dominates the cross-section, which abruptly increases. At each of the resonances, σv depends only marginally on the DM mass m S and, therefore, we have an approximately constant thermally-averaged maximal cross-section (a small m S -dependence arises only at very large values of m S ). This contribution was studied in detail in ref. [21], where it was shown that the resonant enhancement of the cross-section for m S ∼ m Gn /2 was not enough to achieve the value of σ FO v that gives the correct relic abundance, once values of Λ compatible with LHC exclusion bounds were taken into account.
On the other hand, for m S ≥ m G 1 DM annihilation into on-shell KK-gravitons becomes possible. Depending on the DM particle mass, production of several KK-graviton modes is allowed. This is represented in figure 1 by dashed or dot-dashed lines, where we show the contribution to the DM annihilation cross-section from the channels SS → . More channels open for larger values of m S that however have not been depicted in figure 1, where we have decided to show just the lowestlying ones for the sake of clarity of the plot. Recall that each of the two KK-gravitons can have any KK-number: in particular, it is not forbidden by any symmetry to have SS → G n G m with n = m, as translational invariance in the 5th-dimension is explicitly broken due to the presence of the IR-and UV-branes and the KK-number is not conserved. As it can be seen in the figure, the contribution of each channel to the total cross-section varies with the DM mass. For example, SS → G 2 G 2 (orange, dot-dashed line) dominates over SS → G 1 G 3 (green, dashed line) in a very small range of m S , whereas the latter takes over for large m S . Notice that, although KK-graviton production was considered in ref. [18], the possibility of producing different KK-graviton modes was overlooked there.
In figure 2 (left panel) we plot the different Feynman diagrams that contribute to DM annihilation into on-shell KK-gravitons. Diagram (c) was not considered previously in the literature (see, e.g., ref. [18]). However, it must be taken into account when computing the production of two real gravitons, as the corresponding amplitude is also proportional to 1/Λ 2 , the same order as the two other diagrams. 2 The corresponding Feynman rule can be obtained by expanding the metric up to second order about the Minkowski space-time: Notice that, if a diagram that should be considered at a given order in 1/Λ when computing a given process is absent, then the gravitational gauge-invariance of the amplitude is not guaranteed and the cross-section computation is built over slippery ground from a theoretical point of view. The impact of diagram (c) is shown in figure 2 (right panel), where we compare the total DM annihilation cross-section through and into KK-gravitons including or not the contribution to the amplitude from this diagram, for a particular choice of m G 1 = 8 TeV and Λ = 10 TeV. The solid orange (blue dashed) line is the total DM annihilation cross-section through and into KK-gravitons with (without) diagram (c), whereas the dot-dashed red (dotted green) line is the DM annihilation cross-section into KK-gravitons with (without) diagram (c). It can be seen that, for this particular choice of m G 1 and Λ, the difference between the two computations can be as large a one order of magnitude for m S ∼ 10 TeV.
In figure 3 we eventually show the total contribution of KK-gravitons to σv , summing virtual KK-graviton exchange and KK-graviton direct production with contributions from JHEP01(2020)161  Figure 3. The thermally-averaged scalar DM annihilation cross-section through virtual KKgraviton exchange and direct production of two KK-gravitons, σ G = σ ve,G + σ GG , as a function of the DM mass m S . In all panels, the solid orange (blue dashed) lines represent the total cross-section including (not including) mixed KK-graviton production and diagram (c) contribution. The latter case corresponds to refs. [18] and [21]. In order to appreciate the difference,  figure 2 is an increase of the cross-section, that can be as large as a factor two for some specific choices of Λ and m G 1 . In all panels, the horizontal red dashed line corresponds to the value of the thermally-averaged cross-section for which the correct relic abundance is achieved, σ FO v = 2.2 × 10 −26 cm 3 /s. As it was reported in ref. [21], σ FO v is not achievable through KK-graviton exchange since, even for values of m S such that s ∼ m 2 Gn , the resonant cross-section is way smaller than the required one. This result is general and can be found for any value of Λ and m Gn , not only for the few examples shown in figure 3. On the other hand, as reported in ref. [18], for larger values of m S , when the two on-shell KK-graviton production channels take over, a cross-section as large as σ FO v is achievable and the correct relic abundance can be then reproduced. With respect to ref. [18], the net effect of mixed KK-gravitons production and of diagram (c) is to lower slightly the value of m S for which σv = σ FO v . In figure 3, the red-shaded area represents the theoretical unitarity bound σv ≥ 1/s, where we can no longer trust the theory outlined in section 2 and higher-order operators should be taken into account. Notice that, even if in figure 3 the "untrustable" region seems to be very near to the value of m S for which the correct relic abundance can be achieved, it is indeed at least one order of magnitude away, as plots are shown in bi-logarithmic scale.

Virtual radion exchange and on-shell radion production
Consider now the case of DM annihilation into SM particles through radion exchange and of direct production of two on-shell radions, The analytic expressions for the two relevant radion channels contributing to σ r can be found in appendix D.2, whereas in appendix C.2 we give the radion partial decay widths. It can be seen that radion decay to fermions is proportional to the fermion mass squared, Γ(r → ψ ψ) ∝ m r m 2 ψ /Λ 2 , whilst radion decay to bosons (either scalar or vector ones) is Γ(r → B B) ∝ m 3 r /Λ 2 . Clearly, for radions with O(TeV) mass bosons decay channels dominate over fermion ones. However, the decay to massive or massless bosons is rather different: the radion decays to photons and gluons at the one-loop level and, therefore, these decay channels are suppressed with respect to decays into massive bosons, which proceed at tree level. In summary, the radion decay width is dominated by r → W W, r → ZZ and r → HH (and r → SS if kinematically possible).
The two contributions to σ r are shown in figure 4, where we plot σ r (green, dashed line) as a function of m S and compare it with σ G (orange, dot-dashed line). The sum of σ r and σ G is represented by the solid (blue) line. The input parameters for these plots are: m G 1 = 3 TeV and m r = 1 TeV; Λ = 5 TeV (left panel) and Λ = 8 TeV (right panel). For these particular choices of m G 1 , only a couple of KK-graviton resonances appear in σ G before two KK-graviton production takes over. Again, the red-shaded area represents the theoretical unitarity bound σv ≥ 1/s, where we can no longer trust the theory outlined in section 2, whilst the red dashed horizontal line is σ FO v . We can see that, generically and differently from the KK-graviton case, the correct relic abundance can be achieved by the resonant virtual radion exchange channel for DM masses around m S ∼ m r /2 1 + O m 2 r /Λ 2 . Since the radion decay width is rather small, for allowed values of Λ and radion masses in the TeV range or below, a significant amount of fine-tuning is needed in order to get the resonant behaviour. In the absence of a theoretical framework to explain the specific required relation between m S and m r , we consider difficult to defend this possibility as an appealing scenario to achieve the observed DM relic abundance. On the other hand, as it was the case for the KK-graviton exchange and production shown in figure 3, the target value of σv can be achieved also in the range of DM masses for which radion and/or KK-graviton production dominate the cross-section. For the specific values of m G 1 , m r and Λ shown in figure 4 this occurs through KK-graviton production. We have found that this channel dominates in most of the allowed parameter space, while the contribution of radion production is dominant only near the untrustable region m G 1 ∼ Λ.
In figure 5 we show the values of Λ for which the correct DM relic abundance is obtained in the (m S , m G 1 ) plane. In the left panel we assume that the extra-dimension length is stabilized without introducing the radion field. We can see that σ FO v can be achieved in a significant part of the parameter space through KK-graviton production. In order to obtain the target relic abundance σ FO v for m S < m G 1 , small values of Λ are needed, usually excluded by LHC data (as we will see in the next section). Eventually, resonant virtual KK-graviton exchange is not enough to achieve σ FO v for m S m G 1 for any value of Λ, as it is depicted by the grey region (in agreement with ref. [21]).
In the right panel we consider, instead, that the extra-dimension length is stabilized using the Goldberger-Wise mechanism and we introduce a radion with mass m r = 100 GeV. In this case, it is always possible to achieve the correct relic abundance: either through reso--14 -JHEP01(2020)161 nant radion exchange for m S ∼ 50 GeV (not shown in the plot), through radion production in the region m S ≤ m G 1 or, for m S > m G 1 , through KK-graviton production.

Experimental bounds and theoretical constraints
As we have seen in figure 5, in principle the target relic abundance can be achieved in a vast region of the (m S , m G 1 ) parameter space, for Λ ranging from 10 −1 TeV to 10 5 TeV. However, experimental searches for resonances strongly constrain m G 1 and Λ. We will summarize here the relevant experimental bounds and see how only a relatively small region of the parameter space is indeed allowed.

LHC bounds
The strongest constraints are given by the resonance searches at LHC. In our model we have considered two types of particles that could be resonantly produced at the LHC, the KK-gravitons and the radion. In order to quantify the impact of LHC data in our parameter space, first of all we need to compute their production cross-section at the LHC.
The n-th KK-graviton production cross-section at LHC is given by [44]: In our calculations we use the Parton Distribution Functions (PDF's) f i (x) at Q 2 = m 2 Gn obtained from MSTW2008 at leading-order [45]. Regarding the radion, since theq q r vertex is proportional to the corresponding quark mass, the production cross-section in p p collisions at the LHC is dominated by gluon fusion. The gluon-radion interaction is similar to the gluon-Higgs interaction in the SM. We therefore may use the well-known results obtained for the SM Higgs production [46] rescaling the Lagrangian by a factor 3vC 3 /(2 √ 6Λ), where v is the standard model VEV. The final expression is given by: In figure 6 we show the production cross-sections for Λ = 5 TeV at √ s = 13 TeV, where the solid (orange) line stands for p p → G 1 and the dashed (purple) line for p p → r. It is straightforward to obtain the production cross-sections for a different value of Λ by rescaling this plot. As we can see, the radion production is smaller than graviton production by some orders of magnitude. For this reason, the LHC constraints on the Randall-Sundrum model are dominated by (resonant) KK-graviton searches. The KK-graviton decay channels that provide the stringest bounds on m G 1 and Λ are G 1 → γγ [47] and G 1 → [48]. In figure 7 we plot the functional dependence over Λ and  m G 1 of the cross-section p p → , with σ × BR(G 1 → ) ranging from 10 2 fb (bottom line) to 10 −3 fb (top line). Comparing the theoretical expectation with the experimental bounds on σ(p p → ) it is possible to draw exclusion regions in the (m G 1 , Λ) plane, given by the darker (blue) shaded area. The same can be done using the channel p p → γ γ, represented by the lighter (light red) shaded area. We can see that the stringest bounds on Λ are set by p p → G 1 → γγ. Notice that experimental exclusion bounds are given for m G 1 ≥ 200 GeV, approximately.
In figure 8 we show the statistical uncertainties on the experimental bound on σ(p p → ) (left panel) and σ(p p → γ γ) (right panel), where the yellow and green bands are the bounds at 1σ and 2σ in the (m G 1 , Λ) plane, respectively. It can be seen that for low KK-graviton mass the bounds on Λ suffer from a large indetermination: in this range we can only say that Λ should be larger than some value ranging from 50 to 100 TeV, approximately.

Direct and indirect Dark Matter detection
The total cross-section for spin-independent elastic scattering between dark matter and nuclei reads [24]: where m p is the proton mass, while Z and A are the number of protons and the atomic number. The nucleon form factors are given by Due to the fact that in the excluded region the dominant channel to achieve σ FO v is KKgraviton virtual exchange, the exclusion bounds will show a characteristic striped pattern (as it will be shown in figure 10). We have found, however, that constraints from DD experiments are always much weaker than those obtained at the LHC.
Regarding DM indirect searches, there are several experiments looking for astrophysical signals: for instance, the Fermi-LAT collaboration has analyzed the gamma ray flux arriving at the Earth from Dwarf spheroidal galaxies [51] and the galactic center [52,53],

JHEP01(2020)161
while AMS-02 has reported data about the positrons [54] and antiprotons [55] coming from the center of the galaxy. These results are relevant for DM models that generate a continuum spectra of different SM particles, such as the RS scenario we are considering. Recall that DM annihilation into a pair of SM particles via KK-graviton exchange is dwave-suppressed and, therefore, only the annihilation channels into either KK-gravitons or radions lead to observable signals. Both of them will then decay into SM particles leading to a continuum spectrum. 3 However, current data from indirect detection experiments allows to constrain DM masses below ∼ 100 GeV (provided the annihilation cross-section is not velocity suppressed), while for our case of heavy DM (∼ 1 TeV) the limits on the crosssection are well above the required value σ FO v . Thus, indirect searches have no impact on the viable parameter space (see however ref. [19] for other DM scenarios based on RS).

Theoretical constraints
Besides the experimental limits, there are mainly two theoretical concerns about the validity of our calculations which affect part of the (m S , m G 1 , Λ) parameter space. The first one is related to the fact that we are performing just a tree-level computation of the relevant DM annihilation cross-sections, and we should worry about unitarity issues. In particular, the t-channel annihilation cross-section into a pair of KK-gravitons, σ(SS → G n G m ), diverges as m 8 S /(m 4 Gn m 4 Gm ) in the non-relativistic limit s m 2 S , so it is important to check that the effective theory is still unitary. We estimate the unitarity bound as σ < 1/s 1/m 2 S , showing as a green-meshed area in figure 10 the region in which such bound is not satisfied and therefore our calculation is not fully reliable.
The second theoretical issue refers to the consistency of the effective theory framework: in the Randall-Sundrum scenario, at energies somewhat larger than Λ the KK-gravitons are strongly coupled and the five-dimensional field theory from which we start is no longer valid. We therefore impose that at least m G 1 < Λ to trust our results. 4 Notice that this constraint is general for any effective field theory: since we are including the first KKgravitons in the low energy spectra, for the effective theory to make sense the cut-off scale Λ should be larger than the masses of such states.

Achieving the DM relic abundance in RS
We show in this section the allowed parameter space for which the target value of σv needed to achieve the correct DM relic abundance in the freeze-out scenario, ( σ FO v = 2.2 × 10 −26 cm 3 /s) can be obtained, taking into account both the experimental bounds and the theoretical constraints outlined in section 4.
Our final results are shown in figure 10, where we draw the allowed regions of the (m S , m G 1 ) plane for which σv = σ FO v . In the left panel, we are agnostic about the extra-dimension length stabilization mechanism, and assume that neither the unspecified mechanism nor the radion have an impact on the DM phenomenology, as would be the case JHEP01(2020)161 the extra-dimension length is stabilized with the Goldberger-Wise mechanism, with radion mass m r = 100 GeV. In both panels, the grey area represents the part of the parameter space where it is impossible to achieve the correct relic abundance; the red-meshed area is the region for which the low-energy RS effective theory is untrustable, as Λ < m G1 ; the wiggled red area in the lower left corner is the region excluded by DD experiments; the blue area is excluded by resonant KK-graviton searches at the LHC with 36 fb −1 at √ s = 13 TeV; the dotted blue lines represent the expected LHC exclusion bounds at the end of the Run III (with ∼ 300 fb −1 ) and at the HL-LHC (with ∼ 3000 fb −1 ); eventually, the green-meshed area on the right is the region where the theoretical unitarity constraints are not fulfilled. In the left panel, the allowed region is represented by the white area, for which σ FO v is obtained through on-shell KK-graviton production. In the right panel, in addition to the white area, within the tiny orange region σ FO v is obtained through onshell radion production and virtual radion exchange. The dashed lines depicted in the white region represent the values of Λ needed to obtain the correct relic abundance (as in figure 5 of section 3).
for instance if all the new particles in this sector are heavier than the TeV scale; in the right panel, we take into account the radion and consider the Goldberger-Wise mechanism to stabilize the extra-dimension length. The radion mass in this case can be somewhat smaller than the TeV scale (see section 2.2), and therefore it can be relevant for DM annihilation, as we will discuss below. We show our findings for m r = 100 GeV, but other values of m r lead to similar results. As a guidance, the dashed lines taken from figure 5 represent the values of Λ needed to achieve the relic abundance in a particular point of the (m S , m G 1 ) plane. The color legend for the two plots is given in the figure caption.

KK-graviton contributions
Let's consider first the case in which the relic abundance is obtained through virtual KKgraviton exchange and/or on-shell KK-graviton production (left panel). We can distinguish -20 -

JHEP01(2020)161
two regions of the parameter space: 1. m G 1 > m S In this regime the DM annihilates via KK-graviton exchange to SM particles, only. As we have seen in figure 1, the annihilation cross-section is rather small. The grey shaded area in the plot represents the region of the (m S , m G 1 ) plane for which it is not possible to get σ FO v . Below this region, in principle we could find a value of Λ low enough to reach the target relic abundance via resonant KK-graviton exchange. This is, however, in conflict with exclusion bounds in the (m G 1 , Λ) plane from LHC (see figure 7), represented by the darkest (blue) shaded area In addition to the stringent LHC Run II bounds, if the Λ needed to achieve σ FO v for a given m S is smaller than m G 1 , we can no longer trust the RS model as a viable effective low-energy formulation of gravity (diagonal red-meshed area). Therefore, due to the combination of experimental bounds and theoretical constraints, for m G 1 > m S is not possible to obtain σ FO v , as it was indeed found in ref. [21].
2. m G 1 < m S In this case, although the S S → SM SM channel is still open, the target cross-section is achievable through production of on-shell KK-gravitons, S S → G n G m . Due to the LHC Run II bounds, the region of the (m G 1 , Λ) plane for which we can obtain σ FO v corresponds mainly to the region for which (m G 1 /m S ) 2 1. In this region, the value of Λ needed to reach the freeze-out relic abundance is in the range Λ ∈ [10, 10 4 ] GeV, in agreement with the stringent LHC Run II bounds on Λ for relatively low m G 1 . At large values of m S the theoretical unitarity bound discussed in section 4.3 is relevant and, therefore, m S cannot be much larger than 10 TeV (vertical green-meshed area). Eventually, the white area represents the region of the parameter space for which the freeze-out scenario can produce the correct DM relic abundance. Notice that most of this region could be tested either by the LHC Run III 5 (with expected 300 fb −1 ) or by the High-Luminosity LHC (with nominal 3000 fb −1 ), as shown by the dotted lines depicted in the figure. Typical values for m S , m G 1 and Λ in the region that would still be allowed after HL-LHC are m S ∈ [3,15] TeV, m G 1 < 1 TeV and Λ > 10 3 TeV (although a tiny region around m S ∼ 10 TeV with m G 1 as large as few TeV with Λ ∈ [10, 100] TeV could also be viable).
The wiggled dark shaded (red) region in the lower left corner is the bound imposed by XENON1T. The peculiar shape of the bound is a consequence of the resonances in the DM annihilation channels via virtual graviton exchange (see figure 9). We can see that the DD bounds are much weaker than those from the LHC.

Radion contribution
Let's consider now the case in which, in addition to virtual KK-graviton exchange and/or on-shell KK-gravitons production, DM could also produce virtual or real radions (right JHEP01(2020)161 panel). To make easy the comparison with the previous situation, we again consider two regimes: It is always possible to achieve the correct relic abundance through resonant virtual radion exchange and on-shell radion production (see figure 4). In the right plot of figure 10 the former would occur for m S = 50 GeV, outside the range depicted in the figure. Being the radion width extremely narrow, this is possible only in presence of a significant fine-tuning of the DM mass m S and of the radion mass, 2m S ∼ m r . In the absence of a theoretical motivation for such a relation between two, in principle, uncorrelated parameters, we consider this mechanism to achieve the target relic abundance not natural. In the region considered in the plot, the relic abundance can be also achieved through production of on-shell radions for very low values of Λ. This region is represented by the orange (lightest) shaded area. Most of this region, however, is excluded when asking Λ to be larger than m G 1 , as one can see by the diagonal red-meshed area in the plot, Λ < m G 1 . After taking into account the LHC Run II bounds and the limit of validity of the RS model as an effective low-energy theory, a tiny orange-shaded region at m S ∼ 4 TeV, m G 1 ∼ 5 TeV and Λ ∈ [5, 10] TeV is still not excluded. Most of it will be tested with the LHC Run III.
2. m G 1 < m S Since the real KK-graviton production channel, once kinematically open grows very fast as (m S /m G 1 ) 8 (see figure 4), it easily dominates the cross-section. Therefore, in this region of the parameter space there are no significant differences with respect to the case in which the radion is absent, discussed in section 5.1.

Remarks about other setups
In this paper we have focused on the original RS model, in which all the SM particles (and also the DM in our case) are localized on the IR-brane. In the absence of graviton brane localized kinetic terms (BLKT's), within this setup all the SM and DM fields couple to the full tower of KK-graviton excitations with universal strength, Λ −1 . As we have seen, the strong bounds from LHC Run II lead to quite large allowed values of Λ ( 10 TeV), which somehow reintroduce a little hierarchy problem. However many other different configurations have been studied, allowing for some of (or all) the SM fields to propagate in the bulk; for instance, placing gauge bosons and fermions in the bulk has the potential to also explain the hierarchy of fermion masses. Moreover, these extradimensional scenarios can be interpreted as strongly-coupled models in four dimensions (see ref. [18] for details of this duality). Several of the above possibilities have been already analyzed in the context of gravitymediated DM that we are addressing, including DM candidates of various spins (0,1/2 and 1). The idea is that the propagation of SM fields in the bulk and the introduction of BLKT's can reduce suitably the coupling of the SM particles to the KK-gravitons, relaxing the LHC bounds and allowing for lower values of Λ which would then satisfy the original motivation of RS models for solving the hierarchy problem. Although to study in detail -22 -JHEP01(2020)161 these alternative RS scenarios is beyond the scope of this paper, we want to comment in this section about the impact of our results on such other models.
In ref. [21], besides the scenario considered here with all SM and DM fields localized in the IR-brane, two additional benchmark models were studied: 1) SM gauge bosons in the bulk with third generation quarks confined in the IR brane, and all other SM fermions localized close to the UV-brane, so that their couplings to the KK-graviton modes are negligible, and 2) SM fermions localized at various places in the bulk to explain the observed fermion masses and SM gauge bosons propagating also in the bulk. In all scenarios, the Higgs field should remain close to the IR-brane to solve the hierarchy problem, and the DM is also assumed to be localized on the IR-brane. While in none of these setups it was possible to obtain the correct relic density for scalar DM through virtual KK-graviton exchange, the authors did not consider the annihilation channel SS → G n G m nor SS → rr. Since these channels will occur with the same cross-section as in the IR-brane model we analyzed in this paper, it is clear that also in the cases considered in ref. [21] it would be possible to get the target value σ FO v when m S > m G 1 . Actually, it would be easier than in the case considered here, as the LHC bounds on Λ are weaker.
In ref. [19] two additional setups where analyzed and also confronted with indirect bounds from astrophysical data: model A, which addresses the hierarchy problem with the Higgs and DM localized on the IR-brane and the SM matter on the UV-brane, and model B (that gives up the hierarchy problem) where only DM is localized on the IR-brane while the SM matter and Higgs fields are confined to the UV-brane. In both cases, SM gauge bosons propagate in the bulk, so that there is a hierarchy of couplings of the KK-graviton modes, being of order Λ −1 for DM (and the Higgs field in model A) but conveniently suppressed for gauge bosons and negligible for SM matter fields (and the Higgs in model B). As a consequence, the standard radion and KK-graviton searches at LHC do not apply to these models and other searches should be re-interpreted to obtain bounds. Therefore, much lower values of Λ and m G 1 would still be allowed and it should be possible to achieve the correct relic abundance for DM masses in a wider range, from few GeV to TeV, in agreement with our results in figure 5.
In the dual picture of the RS model, the radion is dual to the dilaton, the Goldstone boson from dilatation symmetry in 4D. The dilaton couplings are fixed by scale invariance, and turn out to have the same structure as the radion couplings at linear order. In refs. [34,56], the case in which DM couples to the SM only through a dilaton was studied The authors found that the correct relic abundance can be achieved for light dilaton and DM, since collider bounds from dilaton searches are weaker than for the KK-graviton modes (the dilaton production cross-section is about two-three orders of magnitude smaller than the KK-graviton one, as we can see in figure 6). However, as we are studying a consistent gravitational theory and not only the SM plus a dilaton field, the much stringent bounds from KK-gravitons searches do apply.

Conclusions
In this paper we have explored the possibility that the observed Dark Matter component in the Universe is represented by some new scalar particle with a mass in the TeV range.

JHEP01(2020)161
This particle interacts with the SM particles only gravitationally (in agreement with nonobservation of DM signals at both direct and indirect detection DM experiments). Although this hypothesis would, in principle, mean that the interaction with SM particles is too feeble to reproduce the observed DM relic abundance, we show that this is not the case once this setup is embedded in a warped extra-dimensional space-time, along the ideas of the Randall-Sundrum proposal of ref. [15]. We consider, therefore, two 4-dimensional branes in a 5-dimensional AdS 5 space-time at a separation r c , very small compared with present bounds on deviations from Newton's law. On one of the branes, the so-called "IR-brane", both the SM particles and a scalar DM particle are confined, with no particle allowed to escape from the branes to explore the bulk. In this particular extra-dimensional setup, gravitational interaction between particles on the IR-brane, in our case between a scalar DM particle and any of the SM particles, occurs with an amplitude proportional to 1/M 2 P when the two particles exchange a graviton zero-mode, but with a suppression factor 1/Λ 2 when they do interact exchanging higher KK-graviton modes. Since Λ can be as low as a few TeV (due to the warping effect induced by the curvature of the space-time along the brane separation), clearly a huge enhancement of the cross-section is possible with respect to standard linearized General Relativity.
Using this mechanism, it was studied in the literature if the observed relic abundance in the Universe can be obtained through resonant KK-graviton exchange via σ(DM DM → G n → SM SM) (for any spin of the DM particle), showing that taking into account the LHC bounds on Λ as a function of the mass of the first KK-graviton, m G 1 , it is impossible to achieve the target value of the thermally-averaged cross-section σ FO v for any value of m DM if the DM particle has spin 0 or 1/2 [21]. In refs. [18][19][20]24] it was however shown that, for DM masses larger then the KK-graviton mass, another annihilation channel opens, namely DM annihilation into two (identical) KK-gravitons, σ(DM DM → G n G n ). In this paper, we have studied the possibility that this channel may give a cross-section large enough to attain the observed relic abundance, for the particular case of a scalar DM particle with mass m S . We have indeed found that this is the case and that the region of the parameter space for which σ v ∼ σ FO v is typically at m S of the order of a few TeV, compatible with present direct production searches at the LHC. In the references above some effects were overlooked, though. In particular, a quadratic interaction of the DM particles with KK-gravitons (i.e. the existence of a S S G n G m vertex when expanding the metric up to second order about the Minkowski metric) was not considered. This amplitude is of the same order in 1/Λ as the t-and u-channel contributions to σ(DM DM → G n G n ) considered in the literature and, by increasing the cross-section at large value of the DM mass, lowers the value of m S needed to achieve the relic abundance at fixed value of m G 1 . The same effect is also induced by the possibility of the DM particles annihilating into different KK-gravitons, σ(DM DM → G m G n ), something allowed since translational invariance along the 5-th dimension is explicitly broken by the presence of the branes. This was also overlooked in the existing literature. These effects and their impact have been discussed extensively in section 3 and appendix D.
After having computed the relevant contributions to the cross-section, we have scanned the parameter space of the model (represented by m S , m G 1 and Λ), looking for regions -24 -

JHEP01(2020)161
in which the observed relic abundance can be achieved. This region has been eventually compared with experimental bounds from resonant searches at the LHC Run II and from direct and indirect DM detection searches, finding which portion of the allowed parameter space is excluded by data. Eventually, we have studied the theoretical unitarity bounds on the mass of the DM particle and on the validity of the RS model as a consistent lowenergy effective theory. Our main result is that a significant portion of the (m S , m G 1 ) plane where m S > m G 1 can reproduce the observed relic abundance, for values of Λ ranging from a few to thousands of TeV and m S ∈ [1, 10] TeV. Unitarity bounds put a (theoretical) upper limit on the mass of the DM particle and, interestingly enough, most part of the allowed parameter space could therefore be tested by the LHC Run III and by the proposed High-Luminosity LHC.
In the presence of a Goldberger-Wise mechanism to stabilize the separation between the two branes, the radion r is expected to be light, m r O(TeV), and DM can also annihilate into SM particles via the exchange of a virtual radion and, for m S > m r , two DM particles can also produce directly two on-shell radions. This has been studied in detail in section 3.2 and appendix D.2. Since, contrary to the KK-graviton mass (strongly related to Λ in the RS setup), the radion mass is in practice a free parameter of the model (depending on the unknown details of the scalar potential in the bulk and of some brane-localized terms), it is possible to achieve σ FO v for any value of m S and m G 1 , even in the case m G 1 > m S , through the resonant radion exchange channel (at the price of introducing a significant, theoretically unappealing, fine-tuning of the DM mass with respect to the radion mass, 2m S ∼ m r ) or through on-shell radion production. The region for m G 1 > m S , however, is mostly excluded due to the fact that the value of Λ needed to reach the target relic abundance is Λ < m G 1 , a condition that makes untrustable the RS model as a valid effective low-energy theory. Apart from a tiny region for which the two radion on-shell production channel dominates in the cross-section, the rest of the allowed parameter space is similar to that found in the absence of a radion. FPA2014-57816-P, FPA2017-85985-P and SEV-2014-0398.

A Spin 2 massive graviton
The propagator of the n-th KK-graviton mode, with mass m n , decay width Γ n and 4momentum k in the unitary gauge is: where P µναβ is the sum of the polarization tensors s µν (k) (being s the spin): and The tensor P µναβ must satisfy several conditions for an on-shell graviton G µν , in order to reduce the number of degrees-of-freedom to the physical ones:

B Feynman rules
We summarize in this appendix the different Feynman rules corresponding to the couplings of scalar DM particles and of SM particles with KK-gravitons and radions.

B.1 Graviton Feynman rules
The vertex that involves one KK-graviton (with n = 0) and two scalars of mass m S is given by: This expression can be used for the coupling of both scalar DM and the SM Higgs boson to KK-gravitons. The Feynman rule corresponding to the interaction of two SM Dirac fermions of mass m ψ with one KK-graviton is given by: The interaction between two SM gauge bosons of mass m A and one KK-graviton is given by: and Eventually, the interaction between two scalar DM particles and two KK-gravitons (coming from a second order expansion of the metric g µν about the Minkowski metric η µν ) is given by: The Feynman rules for the n = 0 KK-graviton can be obtained by the previous ones by replacing Λ with M P . We do not give here the triple KK-graviton vertex, as it is irrelevant for the phenomenological applications of this paper. The same occurs for the vertices between one KK-graviton and two radions and two KK-gravitons and one radion.

B.2 Radion Feynman rules
The radion field r couples with both the SM and the DM particles with the trace of the energy-momentum tensor, T = g µν T µν . The only exception are photons and gluons that, -27 -JHEP01(2020)161 being massless, do not contribute to T at tree-level. However, effective couplings of these fields to the radion are generated through quarks and W loops, and the trace anomaly.
The interaction between one radion and two scalar fields (either the DM or the SM Higgs boson) is given by: The vertex that involves the radion and two SM Dirac fermions takes the form: and, as in the case of the graviton-fermion-fermion vertex, we have: The interaction between two massive SM gauge bosons and one radion is given by: The Feynman rule corresponding to the interaction between two massless SM gauge bosons and one radion is:

JHEP01(2020)161
where α i = α EM , α s for the case of the photons or gluons, respectively, and (B.14) with x q = 4m q /m r and x W = 4m w /m r . The explicit form of F 1/2 and the values of the one-loop β-function coefficients b are given by [34]: UV = −11 + 2n/3, where n is the number of quarks whose mass is smaller than m r /2.
Eventually, the interaction Lagrangian between the DM and the radion up to second order is given by:

C Decay widths
In this appendix we compute the decay widths of KK-gravitons and of the radion, using the Feynman rules given in appendix B.

C.1 KK-graviton decay widths
The KK-graviton can decay into scalar particles (including the Higgs boson, the DM particle, if the mass of the considered KK-graviton is sufficiently large, and the radion), SM fermions, SM gauge bosons and lighter KK-gravitons. Decay widths of KK-gravitons into SM particles, Γ(G n → SM SM), are all proportional to 1/Λ 2 . In particular, the decay width into SM Higgs bosons is given by: where m n is the mass of the n-th KK-graviton (in the main text, this was called m Gn , but we prefer here a shorter notation to increase readability of the formulae). If m n > 2m S , the n-th KK-graviton can decay into two DM particles: On the other hand, the decay widths of KK-gravitons with KK-number n into lighter KK-gravitons are proportional to 1/Λ 6 , as the triple graviton vertex comes from the third order expansion of the metric about the Minskowski spacetime. For this reason, we have not considered these decays when computing the total KK-graviton decay widths. The same happens for the radion: the coupling of the radion with the gravitons arises from the mixing of the radion with the graviscalar h 55 , that eventually couples with KK-gravitons again with a triple graviton vertex, proportional to 1/Λ 3 . Also in this case the decay width Γ(G n → r r) is proportional to 1/Λ 6 and, therefore, negligible.

C.2 Radion decay widths
The decay width of the radion into scalar particles, either the SM Higgs boson or the DM particle if the radion is sufficiently heavy, is given by: The radion decay width into SM massive gauge bosons reads: whereas the decay width into SM massless gauge bosons is: (C.9)

D Annihilation DM cross section
Since in the freeze-out scenario, DM annihilation occurs at small relative velocity of the two DM particles, it is useful to approximate the Mandelstam variable s as: s ≈ m 2 dm (4 + v 2 rel ) . (D.1) Within this approximation, the different scalar products for processes in which two DM particles S's annihilate into two SM particles X's, with incoming and outcoming momenta S(k 1 ) S(k 2 ) → X(k 3 ) X(k 4 ), become: where k 1 · k 1 = k 2 · k 2 = m 2 S , k 3 · k 3 = k 4 · k 4 = m 2 X .

(D.3)
We will always write the annihilation cross-sections at leading order in this velocity expansion.

D.1 Annihilation through and into gravitons
The annihilation of DM particles into SM particles through virtual KK-graviton exchange occurs in d-wave. In the following expressions, S KK stands for the sum over all KK states: where m n is the mass of the n-th KK-graviton.
The annihilation cross-section into two SM Higgs bosons reads: The annihilation cross-section into two SM massive gauge bosons is given by: As it was shown in ref. [18], for DM particle masses larger than the mass of a given KKgraviton mode DM particles may annihilate into two KK-gravitons. In the small velocity approximation, the corresponding cross-section is: