I. INTRODUCTION

The particle nature of dark matter (DM) is one of the major outstanding problems in particle physics and cosmology. There are a plethora of models that are consistent with observations and well motivated from a particle physics point of view [1]. However, when viewed from an effective field theory perspective, these models leave many questions unanswered such as those related to the electroweak and/or cosmological hierarchy problems.

Recently, it has been pointed out [2] that the swampland program (see [3] for a review) can offer a new take on the cosmological hierarchy problem, leading to a unification of dark energy and DM. In the context of the swampland program, it is natural for there to be a tower of particles with mass close to the energy scale set by the cosmological constant:

in Planck units and with α ∼ Oð1Þ. This is essentially the statement of the (anti-)de Sitter distance conjecture [4].

In [2] it was shown that, when this line of thought is followed to its logical conclusion, one is led to an essentially unique model of our Universe, dubbed the dark dimension. In more detail, it was argued that theoretical and experimental constraints select α ¼ 1=4 and allow for only a single tower of Kaluza-Klein (KK) states. This has led to the idea of viewing DM as excitations in this tower [5]. The light matter consists of the Standard Model (SM) fields in addition to a tower of massive KK gravitons. Each particle in the tower is coupled only gravitationally to the SM and to other particles in the tower. In the presence of the KK graviton tower and the (necessary) gravitational couplings, it was also shown in [5] that a population of KK gravitons will be produced in the early Universe, leading to a natural candidate for DM, which can constitute all of the DM. This led to a concrete unification of DM and dark energy. Towers of particles have also been considered previously in the context of DM models where they have been motivated by their genericity in string theory (see for example the dynamical DM scenario of [6,7]). The dark dimension is then a well-motivated model from a theoretical perspective given its consistency with swampland principles and its potential to provide a new take on the cosmological hierarchy problem and its unification with DM. It also has phenomenologically interesting signatures as well as possible alternative scenarios that were investigated in various recent papers [8][9][10][11][12].

The purpose of this work is to explore potentially observable cosmological aspects of the dark-dimension model. In particular, we focus on studying perturbations in the linear regime within a simplified version of the full model that nonetheless captures the important physics. Remarkably, we find that the model, although borne out of very general string-theoretic principles, is consistent with current observations and provides signatures that are within reach of near-future experiments.

This paper is organized as follows. In Sec. II we review the relevant phenomenology of the dark dimension. In this section we introduce some details of the "kick" velocity of the DM and distribution functions that we use in our analysis. In Sec. III we present our results and discuss the effects that play an important role in the constraint that we get. Then, in Sec. IV we discuss the parameters of the darkdimension model in light of the new constraints derived from cosmology as well as other experimental data. Finally, we give some concluding remarks and discuss some open questions in Sec. V.

II. KK GRAVITONS AS DARK MATTER

Let us begin by reviewing the essential features of the dark-dimension scenario. For more details, we refer the reader to [5]. In this scenario, the Universe has a tower of massive gravitons with spacing m KK ∼ Λ 1=4 ∼ 10 -100 meV:

ð2:1Þ

These are the KK modes of the graviton in a mesoscopic dimension of length l ≲ 30 μm. The particles of the SM live on a brane localized in the extra dimension and their coupling to the graviton KK modes is determined [up to

Oð1Þ numbers] by the 5D equivalence principle. This allows the KK gravitons to be produced from a hot brane in the early Universe. In particular, if we start with an empty extra dimension (i.e., the KK graviton sector) and an SM brane at temperature T i ∼ 1 GeV, we produce the correct abundance of KK gravitons to account for all of the observed DM. The coupling between each KK graviton and the SM is of gravitational strength. Given this weak coupling, the tower states play an essential role in this production mechanism since it would be difficult to produce a sufficient amount of DM with only gravitational-strength couplings to a few particle states without going close to Planckian temperatures. The large number of particles in the tower allows for the production of enough DM at temperatures much lower than Planckian temperatures. Generically, the extra dimension need not be homogeneous [13] and this inhomogeneity will induce a coupling between the different graviton KK modes. In particular, this coupling will allow heavier KK gravitons to decay to lighter ones, again via gravitational-strength interactions. The large number of graviton modes available for this decay enhances this decay width relative to decays to the SM (see [5] for more details). Thus, when massive gravitons decay, they predominantly decay to lighter gravitons rather than SM particles. 1 We parametrize the energy scale of this inhomogeneity by δ • m KK , where δ ∼ Oð1Þ, which determines the violation of KK-number conservation. More precisely, a parent particle with mass m k can decay to two daughter particles with masses m j and

However, in practice we simply take ϵ ≈ δ • m KK . Given our lack of knowledge of the precise spacing of the tower (which will depend on the inhomogeneity), we in general also allow for noninteger values of δ. Finally, while δ captures the wavelength of the inhomogeneity of the extra dimension, we introduce another effective parameter β that captures the amplitude of these inhomogeneities. This in particular controls the rate of decay of KK gravitons within the tower.

Most dark gravitons are produced at the initial temperature when T ∼ T i with a mass m ∼ T i , and it is a very good approximation to ignore further graviton production at lower temperatures. By then, the most important effects are the decays. Our interest here is in decays within the KK graviton tower. These cause the DM occupation numbers in each mode to vary in the course of cosmological evolution and effectively lead to a time-dependent DM mass (see Fig. 1). As the particles decay, the average DM mass decreases with time and the distribution creeps to lower masses while keeping its shape. At a given time t, the peak of the DM mass distribution is at

where the parametric dependence can be deduced from the total decay width [5] and the prefactor is fixed from numerical simulations.

FIG. 1. The mass distribution of KK gravitons at different times in cosmic history. The horizontal axis should be thought of as a discretum with spacing m KK . The distribution starts at heavier masses and moves to lower masses as the gravitons decay to lighter gravitons.

1

The decay to the SM of course also occurs and has been investigated in [14].

The intratower graviton decays also lead to another important effect: DM velocities. When a parent KK graviton decays, it imparts a kick velocity to the daughter particles that can affect cosmological structure formation. Intuitively, these particles can escape gravitational potential wells and lead to less structure in the late Universe. The kick velocity can be easily determined from the kinematics of the problem and using the knowledge that the coupling between KK particles is gravitational. This latter fact implies that the most likely decay products are the most massive ones. In other words, the gravitational nature of the decay means that when a KK graviton of mass m decays to two daughter particles, the most likely outcome is two particles of masses close to m=2 (up to the mass violation ϵ ¼ δ • m KK discussed above). For simplicity, we assume this to always be the case and neglect subdominant effects from asymmetric decays where one of the daughter particles is much heavier than the other. Finally, we work in the nonrelativistic regime since relativistic DM is strongly ruled out and this approximation will be shown to be self-consistent given the experimental constraint derived below.

Under these assumptions, the kick velocity is determined by the mass of the parent particle alone. We consider the decay of a KK graviton of mass m and treat ϵ ¼ δ • m KK as small compared to m. This is a very good approximation (as we will see later) since ϵ ∼ m KK ∼ meV, while m ∼ 100 keV. The velocity of the daughter particles is deduced by equating their kinetic energy to the difference in masses, i.e., δ • m KK . This gives

where the time dependence is implied by the time dependence of m DM . We take the velocity as a function of time to be given by Eq. (2.3). One might worry that we should incorporate the history of kick velocities. The justification for ignoring this history follows from two effects. The first is that velocity redshifts as the Universe expands, whereby v decreases as a -1 . The second is that the kick velocity increases with time as t 1=7 . These two effects, in addition to the fact that we are working in the nonrelativistic regime, imply that we can ignore the history of particle kicks and focus only on the kick velocity that a particle receives in the last decay. See also Appendix A for more details.

To sum up, we have threefoot_0 parameters in our model relevant to the analysis in this paper: m KK , δ, and β. However, the data we use is only sensitive to the combination that gives the velocity in Eq. (2.3), and that is the single parameter we implement and constrain in the bulk of this paper. Afterwards, we will use these constraints along with others to derive bounds on the dark-dimension model itself. These will be discussed in Sec. IV.

At this point, it is worth emphasizing again that we only make use of theory and data in the linear regime. As such, our analysis does not capture the implications our model might have on the σ 8 tension. In fact, in the linear regime, the significance of the σ 8 tension is much lower (see, for example, Fig. 3) and there is not really a puzzle to resolve. That said, dark gravitons cause a suppression in the amount of fluctuations on small scales and go in the right direction to resolve the σ 8 tension. It would be interesting to study this effect more closely by running N-body simulations that capture the essential phenomenology such as the kick velocities discussed above in the nonlinear regime. However, this is beyond the scope of this work. Finally, we mention that the loss in mass of the dark graviton model is reminiscent of that studied in [15] and may have implications for the Hubble tension as well.

A. Kick velocities

We have already argued that the mean velocity of DM in our model behaves as in Eq. (2.3). There is also a typical spread in velocities which is of the same order of magnitude as the mean velocity. These two features of our model are important for phenomenology and we capture them by assuming the following background distribution of DM momenta: where p 0 ¼ mv today ðt=t today Þ 1=7 and we have added an Oð1Þ constant η to parametrize the variance of the distribution in terms of the mean. In all of our numerical simulations we set η ¼ 0.25, which follows from the behavior of the velocity distribution at late times (see Appendix A for more details). The time-dependent normalization of this function is chosen such that the DM energy density redshifts like a -3 . As shown in [5], this is a very good approximation despite the decay of DM particles and the subsequent loss of the daughter particle kinetic energy as the Universe expands. In fact, the fractional energy loss in one Hubble time is

ð2:5Þ

We will see later that we need v ∼ 10 -4 , implying a negligible loss in DM energy density due to the decays. Said differently, the energy loss is an Oðv 2 Þ effect which is subdominant to the OðvÞ velocity effects that we study in this work.

We study perturbations by taking

where the Boltzmann equation determines the evolution of Δ. Rather than tracking the evolution of the full distribution function (2.4) and its perturbation via a Boltzmann hierarchy, we work in the fluid approximation and deal with moments of f 0 and its perturbation. We truncate the Boltzmann hierarchy at Legendre multipole l max ¼ 2 (see, for example, [16]). For a different approach see, for example, the treatment of [17]. The relevant definitions and equations governing the perturbations are the familiar ones and are reproduced in Appendix B for completeness.

III. ANALYSIS WITH COSMOLOGICAL DATA

The DM we consider has a few novel features compared to particle physics models often considered in the literature (see, e.g., [18]). First, our DM starts with small velocity which then increases over time. This is an unusual behavior from the point of view of cosmology, where velocities and temperatures decrease with the expansion of the Universe. Moreover, since all of the dark gravitons participate in decays, it is all of the DM that acquires larger velocities with time rather than just a fraction. In addition, our model does not have a dark radiation component unlike models of decaying DM which often include dark radiation as well. Finally, on a more detailed level, the one-seventh power in Eq. (2.3) may allow future experiments to distinguish this model from other models of decaying DM (see below for more on this).

The features mentioned above play an important role in the phenomenology. The main effect is on the evolution of density perturbations (and the related gravitational potentials). In general, one has to do a full N-body simulation to capture the nonlinear behavior of a model. In this study, we only focus on the linear regime (see also Sec. III D for more details). In this regime, we study effects on the cosmic microwave background (CMB) and large-scale structure of the Universe. That the CMB is affected is easy to see through the line-of-sight formalism where the gravitational potentials enter explicitly in the calculation of the temperature anisotropies. The effect on the large-scale structure that we study can be traced back to a suppression in the amount of structure on small (but still linear) scales. We probe these effects indirectly via the changes they induce in the CMB and weak lensing. We rely on Planck 2018 data for measurements of the CMB and KiDS-1000 for weak lensing (see also Sec. III B below).

Using Planck data, we get a robust constraint on v today :

CL; PlanckÞ: ð3:1Þ

These constraints are derived from two effects that alter the CMB. First, the intrinsic unlensed CMB is suppressed at high multipole due to the presence of the graviton velocity which suppresses fluctuations. The other important effect is due to alterations of the lensing potential which relies on the intermediate structure between us and the surface of last scattering. In addition, as can be seen from Fig. 2, there is a negative correlation between S 8 and v today , as expected. This is due to the intuitive fact that larger velocities would lead to less structure. The Planck measurement of S 8 (and other cosmological parameters) has smaller uncertainty and, when combined with weak lensing data, Planck fixes these remaining cosmological parameters.

Let us now briefly discuss the constraints when adding weak-lensing measurements. As mentioned previously, we only make use of linear scales. Since this cutoff between linear and nonlinear scales cannot be defined in a precise way, the constraints on v today depend on where we make this cut. More details on this is provided in Sec. III D. That said, the 95% CL constraints are where the definition of Δ NL is given in Sec. III D and the subscript 0 indicates that this is the largest k in the linear regime today. The mode numbers shown for Δ NL ¼ f0.5; 1.0g are consistent with those found in [19]. It is easy to understand where this constraint comes from intuitively. The gravitons' kick velocities allow them to escape potential wells that would otherwise cause them to collapse into bound structures. This suppresses the amount of structure on small (but still linear) scales and alters the lensing (convergence) power spectrum that the data measures. A large change that comes with a large value of the velocity parameter v today is disfavored by the data. All in all, and depending on the cut between linear and nonlinear scales, the two data setsfoot_1 provide complimentary constraints in the S 8 and v today plane, as shown in Fig. 2. While we quote the bound for two values of Δ NL here, we will use the conservative constraint, with Δ NL ¼ 0.5 in the rest of the paper.

Finally, we comment briefly on related constraints from nonlinear scales that have been discussed in the literature [20,21]. These constraints come from comparing the known abundance of low-mass subhalos with results from N-body simulations that explore galactic substructure. In these decaying DM models, one considers a long-lived DM candidate. The lifetime of such a particle is taken to be close to the age of the Universe, τ ∼ t H , and thus early structure formation proceeds in the same way as in the Λ cold dark matter (ΛCDM) model. However, once DM particles begin to decay (typically at low redshifts because of the chosen lifetime) the evolution can deviate from CDM. This deviation depends on the recoil/kick velocity. Intuitively, if the kick velocity is large enough, then it can strip mass off halos and subhalos and make the latter especially susceptible to tidal disruption events. As such, fewer subhalos survive to the present day and one would observe a reduced abundance. The reported bounds [21] depend in general on the recoil/kick velocities [for example, for kick velocities of 20 km=s (40 km=s), a lifetime longer than 18 Gyr (29 Gyr) is required for the DM candidate to not be ruled out]. However, for recoil velocities larger than the maximum circular velocity of the subhalos, the bound becomes independent of the kick velocity and depends only on the DM lifetime. For kick velocities this large, all decays lead to the daughter particle being ejected from the DM halo and all that matters is the number of decays, which is set by the lifetime. In this case, lifetimes shorter than ∼40 Gyr are ruled out [21].

Our model is slightly different from models of a single decaying dark particle since our DM decays throughout cosmic history rather than just at low redshift. In addition, as we previously argued, the kick velocity is a function of time, as is shown in Eq. (2.3), unlike the constant kick velocities often studied in the literature. However, despite these differences one may hope that the bounds from the literature can be applied to our model conservatively, although it would be impossible to draw quantitative conclusions without further analysis.

Given the discussion above, we check that our model is not ruled out by observations of subhalo abundances. The natural parameter range for our model has a kick velocity that is of the same order as the bounds in Eqs. (3.1) and (3.2), which we can take to be v today ∼ 10 -foot_2 ∼ 30 km. These velocities are similar to the recoil velocities studied in [21]. As such, we are in the regime where only the lifetime of the dark gravitons matters (as explained above) and the lifetime bound τ ≳ 40 Gyr is a constraint on our model. We now show that our model automatically satisfies this bound.

The dark graviton decay lifetime 4 is given by [5]

A. Effect of the kick velocity on cosmology

In this section we discuss the effects of including a kick velocity of the form (2.3) on the evolution of cosmic perturbations. We are mainly interested in small scales that are nonetheless in the linear regime. As usual, the larger scales evolve just as in ΛCDM (cf., for example, the case of warm DM [22,23] where only small scales differ from the ΛCDM predictions). The inclusion of the kick velocity washes out inhomogeneities on small scales, and in this section we discuss the dynamics of this process. We work in synchronous gauge but the translation to any other gauge should be straightforward.

We start with the equations governing the evolution of dark graviton perturbations (in synchronous gauge):

where an overdot denotes a derivative with respect to conformal time τ and the subscript f stands for "fluid," as we are treating the dark gravitons in the fluid approximation. In writing the above equations, we have dropped terms proportional to the shear since they only play a subdominant role in the following discussion. That said, these terms can be found in Appendix B and have been taken into account in the full analysis. The expressions for the equation of state w and the sound speed c 2 a (which is obtained following Ref. [17]) are given by

2=7

; ð3:7Þ where p 0 , ρ 0 , and p 0 are the constant pressure, energy density, and pseudopressure, respectively, and are defined in Appendix B. In particular, the above equation shows that for small v today , we have c 2 a ≈ 5w=3 and the adiabatic sound speed c 2 a and equation of state w are of the same order. The CDM limit of the above perturbation equations is obtained by taking v today → 0, which implies w; c 2 a → 0. In addition to the above two equations, we need equations for the metric perturbations ḣ and all other perturbations that source ḣ. These are standard and can be found in many treatments of cosmological perturbation theory so we omit a detailed discussion and simply write the schematic form of the ḣ equation:

where the sum is over all of the species present in the Universe.

We now consider the evolution of the dark graviton perturbations δ f and θ f for a particular comoving wave number k. The wave numbers of interest are superhorizon modes deep in the radiation era and enter the horizon in the radiation era. These modes start with the same adiabatic initial conditions as cold DM:

where δ ðiÞ γ is the initial condition for the photon density perturbation and the vanishing initial velocity perturbation is the usual one imposed on CDM in synchronous gauge. The choice of CDM initial conditions is justified in our case because the kick velocity vanishes in the limit t → 0, so that our graviton DM reduces to CDM at very early times. This is also apparent in Eqs. (3.4) and (3.5), which also reduce to the CDM evolution equations at early times. Since we have the same initial conditions and evolution equations as CDM at early times, the density and velocity perturbations of dark gravitons evolve just like CDM until other non-CDM terms in Eqs. (3.4) and (3.5) become important.

In particular, deep in the radiation era, the gravitational potentials are set by the radiation perturbations (photons and neutrinos) and the dark graviton perturbations simply evolve in these potential wells. While they are still superhorizon, the density perturbations grow as δ f ∼ a 2 , just like CDM. However, unlike CDM, the velocity perturbation θ f , is small but nonzero. 5 While θ f vanishes initially, the second term of Eq. (3.5) is nonzero and causes the θ f solution to deviate from its CDM counterpart. The two terms are at the same order and θ f is given by the parametric expression

This value is still suppressed (by the smallness of w) compared to δ f and unimportant at this stage. Intuitively, for superhorizon modes, the pressure is not large enough to affect the evolution of the fluid. That said, as δ f grows, the velocity perturbation also grows and this will eventually alter the evolution of this wave number so that it deviates from the CDM solution. For the values of v today and wave numbers relevant to our analysis, this happens in the matter era. Let us then continue following our perturbation mode through horizon crossing and into the matter era. A particular wave number becomes subhorizon when H < k. This happens in the radiation era for modes and parameters of interest to us, as we mentioned previously. In this case, the evolution of δ f no longer follows a 2 and becomes logarithmic instead, going like ln a. This continues until the Universe transitions into matter domination, where δ f starts growing like a power law δ f ∼ a again. Throughout this time, the velocity perturbation is growing, as is given by Eq. (3.10). This behavior continues until θ f is large enough to compete with the gravitational potential ḣ term in Eq. (3.4). At this point, the pressure (in this mode) can counteract the effect of gravitational collapse and the mode stops growing and starts oscillating instead. As such, this mode ends up today with a lower amplitude than it would have had in the CDM case, i.e., the linear matter power spectrum PðkÞ is suppressed. Smaller wavelength modes deviate from the CDM behavior earlier and thus end up with more suppressed fluctuations compared to longer wavelengths.

In the following sections, we will be referring to the Newtonian gauge gravitational potentials, so we briefly review the physics of the suppression of fluctuations in this gauge. Recall first the Newtonian gauge metric perturbations Φ and Ψ that appear in the 00 and ii components. We work in the limit where we ignore the shear, implying (as usual) that Φ ¼ Ψ. The analog of Eq. (3.8) is now the algebraic Poisson equation that determines the gravitational potential Φ directly in terms of the matter perturbation variables δ i and θ i . As such, we focus our discussion on the matter perturbations only. In Newtonian gauge, these are governed by the following equations:

where we are ignoring shear contributions as in the previous discussion in synchronous gauge. Like before, we are interested in modes that enter the horizon during the radiation era. As in the case of the synchronous gauge above, the evolution of these modes splits into two regimes (see, for example, [24]). In the first regime, we focus on horizon entry in a radiation-dominated Universe. This is followed by a second regime which tracks subhorizon evolution through matter-radiation equality. We start by reviewing the evolution of these modes in the CDM limit (obtained by taking w; c 2 a → 0) before moving on to the discussion of dark gravitons. In the standard ΛCDM model, the CDM density perturbation remains constant for superhorizon modes in the radiation era. After horizon entry in the radiation era, CDM perturbations grow logarithmically (proportional to ln a) and then linearly (proportional to a) when the Universe becomes matter dominated. Similar to the matter perturbation, the gravitational potential Φ also remains constant for superhorizon modes in the radiation era. It then decays rapidly in the radiation era (proportional to a -2 ) when it becomes subhorizon. Once the Universe becomes matter dominated, the amplitude of this Φ mode remains constant (until the dark energy epoch, which is beyond our interest here).

Let us contrast the above behavior to that of our model where w; c 2 a ≠ 0. In particular, the θ f Eq. (3.12) contains a term proportional to the density perturbation which grows with time. This term increases for two reasons: because the δ f perturbation is growing, and because the kick velocities cause the sound speed to increase. As such, there comes a time (usually during the matter era for the modes of interest) where this term can interfere with the standard evolution of θ f due to Hubble expansion and the gravitational potential. When this happens, it causes the θ f and (in turn) δ f perturbations to oscillate instead of grow. As such, these density perturbations end up with a smaller amplitude today (compared to their CDM counterparts) and the matter power spectrum is suppressed on these scales. Relatedly, the gravitational potentials, which are determined algebraically in terms of δ i and θ i , are also suppressed.

In summary, the growth of the velocity perturbation for graviton DM leads to a suppression in fluctuations on small scales. The details of this mechanism are quantitatively different from other models relying, for example, on the t 1=7 factor in Eq. (2.3). This will distinguish signals of dark gravitons from other mechanisms that also suppress the power spectrum (for example, warm DM) allowing future experiments to potentially differentiate the dark-dimension scenario from other models.

The cosmological effects we discuss in this paper can all be traced back to this suppression in density (and, relatedly, gravitational potential) perturbations on small scales compared to their ΛCDM counterparts. Ultimately, this leads to a suppression in the matter power spectrum and this may help alleviate tensions between large-scale structure and CMB data within the ΛCDM paradigm. That said, we do find evidence of this in the linear regime, and a more indepth analysis into the nonlinear regime is required to shed light on these tensions. We now turn to a brief discussion of these effects.

B. Data sets

In order to compare our model to cosmological data, we modify the publicly available Boltzmann solver CLASS [25] and interface it with the sampler MontePython [26,27]. The modification adds one additional parameter v today that we vary alongside the usual ΛCDM parameters fA s ; n s ; θ Ã ; Ω b h 2 ; Ω c h 2 ; τ reio g. We use the following data sets in various combinations.

(1) Planck 2018: We use the DR3 version of publicly available likelihood code released by the Planck Collaboration. In particular, we use high-and low-l temperature and polarization data as well as the Planck lensing measurements [28][29][30].

(2) Baryon acoustic oscillations (BAO): We use the combined galaxy and Lyman-α BAO likelihood released in [31]. This likelihood uses data from eBOSS DR14 galaxy BAO measurements and autoand cross-correlation of Lyman-α absorption with quasars [32,33]. (3) KiDS-1000 weak lensing measurements: We use the cosmic shear measurements of KiDS-1000 [34] in the form of E-and B-mode band powers. The bandpower spectra are implemented in a MontePython likelihood [35] and can be used for Markov chain Monte Carlo simulations in a KiDS Cosmology Analysis Pipeline environment [36]. The power spectra are obtained by correlating cosmic shear measurements in five tomographic redshift bins for a total of 15 independent cross-or auto-correlations. We modify the publicly available likelihood to make use of only the linear scales. This is discussed further in Sec. III D. We perform a Markov chain Monte Carlo likelihood analysis, where the sampling is done using the Metropolis-Hastings algorithm, except when deriving constraints from the KiDS-1000 likelihood only where we use nested sampling [37][38][39]. Finally, when running chains using only the KiDS-1000 data, we impose the priors listed in Table I that match those used in [35], for example, as well as other KiDS analyses. In addition, we fix the ΛCDM parameters n s and τ to their Planck measured values as KiDS-1000 data does not constrain them.

C. Effects on the CMB

The CMB can be potentially altered in two ways by the dark graviton kick velocities. First, there is a possible change to the intrinsic unlensed CMB. Neglecting largescale effects (such as the late integrated Sachs-Wolfe effect), this will happen on scales that begin to deviate from their ΛCDM behavior at or before recombination (z ≈ 1100). This happens at very small scales on the order of 10-20 Mpc -1 (recall that larger scales deviate from ΛCDM at later times) and in fact much smaller than those accessible to current observations.

The second important effect is gravitational lensing. Here, we can divide the story into two parts: the lensing power spectrum itself and the lensing of the CMB. The lensing potential quantifies the deviation of light geodesics as they approach an observer and can be calculated by a line-of-sight integral:

where τ is conformal time, χ is the comoving angular diameter distance, and Φ is the Newtonian potential. It is clear that the lensing potential will incur large changes when the Newtonian potential Φ varies. In particular, in our case, where the Newtonian potential in some modes decays at late times, we expect to observe a smaller value of ϕ and its variance on those scales. A large suppression in the lensing power spectrum is disfavored by the Planck lensing likelihood and this is partially responsible for our model constraints.

The lensing potential defined in Eq. (3.13) is also responsible for lensing of the CMB (see, for example, [40,41]). Schematically, the observed lensed temperature anisotropy Θð nÞ ≡ ΔTð nÞ=T is related to its unlensed counterpart Θð nÞ by Θð nÞ ¼ Θð n þ ∇ϕÞ ≈ Θð nÞ þ ∇ϕ • ∇ Θð nÞ ð3:14Þ

in the limit of weak lensing. A smaller lensing potential leads to less distortion of the CMB and sharper peaks in the power spectra. For v today ∼ 10 -3 , this effect contributes an Oð1%Þ in the power spectrum. Finally, we discuss how the above effects change the goodness of fit by appealing to the Planck 2018 likelihoods. To do that, we compare two models: one with v today ≪ 10 -3 (the ΛCDM limit) and the other with v today ≈ 10 -3 . This latter one is at the boundary of the 95% CL region shown in Eq. (3.1). We use the full TTTEEE likelihood as well as versions with the l range restricted to l < 800 and l < 1600 to understand what region of the CMB drives these constraints. The results are presented in Table II. We see that a change in v today (while keeping all other parameters fixed) affects the CMB on all scales. In addition, the fit to the reconstructed lensing power spectrum also deteriorates, as shown by the Planck lensing likelihood. Altering the cosmology cannot improve the CMB fit by much and the largest improvement is to the lensing power spectrum. Nuisance parameters, on the other hand, improve the CMB fit for the most part. Altogether, the overall fit remains slightly worse than that with lower

D. Effects on weak lensing

The observed image of a galaxy can be gravitationally lensed by intervening cosmic structure along the line of sight to the galaxy. This can cause two types of distortions of the original image (for a review see [42]). The first is an isotropic stretching of the image, which is quantified by the convergence. The second is an anisotropic stretching altering the shapes of galaxies and is quantified by the shear. This decomposition of the shear field is exactly analogous to the decomposition of the CMB temperature anisotropies into curl-free E-modes and divergence-free B-modes.

In the case of weak lensing, the B-mode autocorrelation spectrum is expected to vanish in the absence of systematics. As such, we focus our discussion on the E-mode (convergence) power spectrum, although we use the Eand B-mode power spectra supplied by the KiDS-1000 Collaboration [34,35] in our analysis. In the Limber approximation, the convergence power spectrum can be expressed in terms of a line-of-sight integral over the matter power spectrum multiplied by two weight functions q μ ðχÞ (see also [43,44] for more details):

where χ H is the comoving distance to the horizon and P m is the matter power spectrum. The window functions quantify the lensing efficiency and contain information about the galaxy redshift distribution. They are given by

where n μ ðχÞ is a normalized source redshift distribution. The presence of two weight functions is because we are computing a two-point correlation function. The indices μ, ν indicate that, depending on the chosen weight functions, one can compute the convergence power spectrum between different source populations. For instance, in the KiDS-1000 analysis that we follow here, the galaxy sample is divided into five tomographic redshift bins and one considers the 15 independent convergence correlations between these bins. In this case, there are five choices for each weight function and they quantify the galaxy redshift distribution in each tomographic bin.

The original KiDS-1000 analysis makes use of nonlinear approximations to the matter power spectrum P m (for example, using Halofit [45,46] or HMCode [47,48]). However, our model is expected to be different from ΛCDM on nonlinear scales and the prescriptions for deducing the nonlinear power spectrum from the linear one (as is done in Halofit or HMCode) need to be modified. We leave this for future work and focus instead on the linear regime approximating P m by the linear power spectrum P lin m . In addition, we need to ensure that we only use information from linear scales and we do this by restricting the multipole range in the convergence (and shear) power spectra according to the following prescription.

We consider a comoving scale k at redshift z to be nonlinear if it is larger than k NL , where the latter is given by the following criterion:

where P lin m ðk; zÞ is the linear power spectrum and we take three values for Δ NL ∈ f0.5; 1.0g. Doing this allows us to quantify the dependence on k NL since the cutoff between linear and nonlinear scales is not sharp. We emphasize that k NL depends on redshift since the linear power spectrum grows with time. In practice, the above restriction to linear scales means that we only use certain low multipoles for each redshift and discard all of the higher ones. When the spectrum in question is a cross-correlation between two redshift bins, then we take the more conservative choice for the multipole range (i.e., only using the data points that are TABLE II. Breakdown of -2 ln L contributions in going from one best-fit model to another along a path where we change one set of parameters at a time. The two best-fit models are (i) a v today ≪ 10 -3 best-fit model (the ΛCDM limit) and (ii) the best-fit model with a fixed v today ≈ 10 -3 . Fitting is done to full Planck data in addition to the BAO data set discussed in Sec. III B. The table is best read from left to right with changing parameters indicated in the headers and the differences in -2 ln L in the table. The first and last columns show the value of -2 ln L. The benchmark value v today ≈ 10 -3 is chosen since it lies within the 95% CL constraint using Planck data (see Fig. 2). linear at both redshifts). We refer to this restricted data set by the name "Linear K1K" and use it instead of the full KiDS-1000 data when constraining our model. Finally, we note that for Δ NL ∈ f0.5; 1.0g, we get k NL ∈ f0.09; 0.17g Mpc -1 at z ¼ 0, which is roughly when the linear approximation begins to breakdown. Constraints on v today are easily understood by considering Eq. (3.15), which clearly shows that a suppression in the amount of structure [i.e., P m ðk; zÞ] leads to a suppression in the convergence power spectrum. We show an example of the convergence power spectrum in the ΛCDM limit (v today → 0) and the benchmark v today ≈ 10 -3 in Fig. 3. Note that C l 's measuring correlations between higher-redshift bins deviate more from ΛCDM. This is due to the fact that the lensing power is a cumulative effect given by a line-of-sight integral [see Eq. (3.15)]. Light from the higher-redshift tomographic bins travels a longer distance and this gives a bigger difference in the cumulative lensing effect between our model and ΛCDM.

Likelihood

IV. MODEL PARAMETER CONSTRAINTS

In this section we set the cosmological constraint on the velocity of dark gravitons today in the context of our model and other experimental bounds. In particular, we discuss what range of model parameters are compatible with current experimental constraints.

Let us begin by summarizing the constraints we have. First, we have the bound on the velocity today derived in this paper:

ðthis work; linear cosmologyÞ: ð4:1Þ

Next, we have the astrophysical bound derived in [14],

Finally, we also have the bound on the lightest particle in the KK tower from fifth-force experiments [49,50]:

ðfifth-force experimentsÞ:

ð4:3Þ

In the last inequality, we included an Oð1Þ number γ that is the ratio of the mass of the lightest KK particle to m KK , which is the spacing of the tower at asymptotically large KK numbers.

FIG. 3. KiDS-1000 data overlayed with the C EE l best-fit power spectrum calculated for ΛCDM and from our model for a fixed value of v today ≈ 10 -3 . The grey data points have been masked in the analysis since they lie in the nonlinear regime. See Sec. III D for more details.

Before discussing the constraints these imply on our model parameters, we briefly review the parameters themselves and their physical meaning.

(1) λ: This is a measure of a suitably average wavefunction overlap between SM fields that live on a (3 þ 1)-brane and KK gravitons that live in an effectively 5D spacetime. This parameter sets the strength of the coupling between the dark gravitons and the SM. As stated earlier, this coupling is of gravitational strength and is determined parametrically by M Pl , but the wave-function overlap changes the value of the coupling by an Oð1Þ coefficient that we call λ.

(2) β: This is a measure of the amplitude of inhomogeneities in the extra dimension. This parameter sets the strength of the coupling between particles in the dark graviton tower. Intratower decays are faster for larger β.

(3) δ: This describes the wave number of the violation of KK momentum in dark graviton decays (as discussed in Sec. II). ( 4) m KK : This is the spacing between KK particles at large KK quantum numbers and, physically, it is the inverse of the size of the extra dimension. ( 5) γ: This is the ratio of the mass of the lightest particle in the tower to m KK , the asymptotic mass spacing. We note that m DM is not an independent parameter in our model as it can be derived from the remaining parameters. Generally, the parametrics of m DM are set by equating the decay rate to the Hubble rate [as shown in Eq. (2.2) and discussed in Sec. II]. For convenience, we repeat that equation below:

We can get a lower bound on m DM today using Eqs. (

:

In summary, our model has five independent parameters ðλ; β; δ; m KK ; γÞ and two derived parameters ðv today ; m DM Þ that are determined in terms of the free parameters. In addition, we have three constraints given by Eqs. (4.1)-(4.3). These rule out a portion of parameter space and the goal is to determine the largest region that is compatible with all constraints.

The strongest inequalities we can derive involve weighted products of the parameters β and δ. Given values for the remaining parameters, these inequalities will lead to constraints in the ðβ; δÞ plane. The first such inequality is easy to find by rearranging Eq. (2.2) and using the lower bounds on m KK from Eq. ( 4.3) and m -1 DM from Eq. (4.2). This gives

Next, we can also find an upper bound on a weighted product of β and δ. This can be obtained by expressing the DM mass in terms of the velocity (as well as δ and β) and using the inequalities from Eqs. (4.1) and (4.4). This gives

The remaining combination of inequalities is given by Eqs. (4.1) and (4.2), but it is easy to check that this leads to a bound on m KK that is much weaker than the experimental constraint (4.

3) [at least for δ; λ ∼ Oð1Þ]. These last two inequalities provide the strongest relations between the dimensionless parameters of our model. They can be combined (e.g., by eliminating one of β or δ) to obtain That said, it is more illuminating to work with the bounds on the weighted products themselves [Eqs. (4.5) and (4.6)], so we proceed to do that.

To get insight into the allowed values for β and δ, we have to choose values for the parameters λ and γ. The natural region for the parameter λ has been discussed in [14]. There, it was found that λ ∼ Oð1Þ is allowed by astrophysical observations and large λ (which is anyhow unnatural from the point of view of our model) is ruled out. As such, in this work, we take two benchmark values for λ ¼ f0.1; 0.5g. For the parameter γ, we take three values f0.5; 1.0; 2.0g. We show the allowed regions in Fig. 4.

Finally, using the above constraints we can get a lower bound on the size of the extra dimension (equivalently, an upper bound on m KK ). This can be easily obtained using the upper bound on v today from Eq. (4.1) and the upper bound on the mass of the dark gravitons from Eq. (4.4). These give In particular, note that the overall inequality reduces to that in (4.7).

V. OUTLOOK AND CONCLUSION

In this work, we studied a model of DM motivated by swampland considerations. In this model, dubbed the "dark dimension," DM is composed of a tower of massive spin-2 particles (gravitons) with spacing on the meV scale. Particles in the tower are populated in the early Universe and decay to lighter particles as the Universe evolves. Each of these decays gives a kick velocity to the daughter particles, which has an effect on cosmology. It is precisely this effect that we investigated in this work.

We argued why the DM kick velocity is time-dependent and used the fluid approximation to study how cosmological perturbations evolve in the presence of this kick velocity. We found, as expected, that the kicks can inhibit the growth of structure on small linear scales and that this leads to a constraint on the velocity today. The constraint originates from effects on the CMB and weak lensing. That said, velocities in the natural range for our model are viable and may be observed in future experiments.

For our analysis, we focused on the linear regime. We leave a study of the dark-dimension model in the nonlinear regime for future work. Exploring this question can lead to better understanding of the model and improved constraints.

Related to the above, large-scale structure surveys will soon present us with a lot more data [51]. It would be interesting to see what that data means in the context of the dark-dimension model. As emphasized previously [2,5,14], the natural parameter range for this model is close to current limits. This view is also reinforced by the analysis in this paper. With new data, the limits will improve and could lead to a detection or ruling out of the model with potential implications for our understanding of quantum gravity.

APPENDIX A: VELOCITY DISTRIBUTION

Let us justify the velocity distribution we use in Eq. (2.4). Since we are interested in the nonrelativistic limit, the problem is reduced to the addition of velocity distributions taking into account the vector nature of the particle velocities. Moreover, we are simplifying our problem by assuming that the kick velocity depends only on time. In this appendix we use another simplification, which is to treat the velocity kicks as happening at discrete time intervals. Each of these time intervals corresponds to a Hubble time where roughly all of our dark graviton particles that we had originally have decayed.

To start with, let us ignore the expansion of the Universe and find the evolution of the velocity distribution in flat space ignoring gravity. Let us assume that we start with a collection of particles, all of which are at rest so that the initial velocity distribution is given by χðv

where we emphasize that all distributions we consider are isotropic and depend only on the magnitude of the velocity vector v ¼ j⃗ vj. In fact, we consider the distributions as functions of v 2 rather than v, which we will see is a more convenient choice. Finally, we point out that the definition of χ is normalized such that R dv 2 χðv 2 Þ ¼ 1. At the first time step (the analogue of a Hubble time in the cosmological scenario), all of these particles decay and the daughter particles all get a nonzero velocity. However, given that the magnitude of the kick velocity depends only on time, all of the daughter particles will have velocities of the same magnitude but different directions. The velocity distribution then becomes

where v 1 is the kick velocity obtained in going from t 0 to t 1 .

At the first time step, we see that the particles are no longer at rest. Instead, their velocity distribution is given by a δ function.

At the next time step, in going from t 1 to t 2 , again all particles decay and this imparts a kick velocity ⃗ v 2 to be added to the original velocity of the parent particle ⃗ v 1 . The magnitude of the resulting vector ⃗ v 3 is given by

where θ 12 is the angle between ⃗ v 1 and ⃗ v 2 . Since the decay is isotropic, the distribution of the cosine term is uniform, cos θ 12 ∼ Uniformð-1; 1Þ; ðA4Þ so that v 2 3 is also distributed uniformly between jv 1v 2 j 2 and jv 1 þ v 2 j 2 :

where Θ is the Heaviside function that vanishes for a negative value of its argument and is unity otherwise. The fact that the expression (A3) gives v 2 is the motivation for taking χ to be a function of v 2 . We see that the distribution becomes wider (going from a δ function to Uniform) as we step forward in time, and this is expected. However, our interest is in the late-time behavior of this distribution. So we have to iterate the above procedure for many time steps and describe the late-time result. We will do this numerically but, in order to do that, we need a formula to obtain χðv 2 ; t þ ΔtÞ from χðv 2 ; tÞ for a general function χ. We do this by relying on this last observation which says that a δ function in velocity becomes a uniform distribution at the next time step. All we have to do is decompose any distribution of interest into a "sum" over δ functions and replace each of these by a uniform distribution.

We now take a general velocity distribution χðv 2 ; tÞ and deduce the time-evolved distribution χðv 2 ; t þ ΔtÞ by summing over contributions from all points on χðv 2 ; tÞ, each of has with an amplitude like that given in Eq. (A5). Let us call the kick velocity v k . We then immediately have

As a quick check, we can test the above formula by taking χð v2 ; tÞ ¼ δðv 2v 2 Ã Þ, and we recover Eq. (A5). Using the above formula, and taking v k to increase as a function of time, it is now straightforward to show (numerically, as we do) that χðv 2 ; tÞ has a peak around v k and a width of order v k as well. This is the late-time behavior of the distribution χðv 2 ; tÞ that we are after.

We can now use the above equation to get an intuition for how the velocity distribution changes in our model. However, before doing that, we need to include the cosmological expansion which also causes a redshift of the velocity distribution between time steps.

In order to include the cosmological expansion, first observe that the effect of redshift is to transform a distribution χ in the following way:

where a 0 ¼ aðt 0 Þ. The overall factor is important to preserve the normalization R χdv 2 ¼ 1 at all times. We can then see that if we were to work with the quantity w ¼ aðtÞv and consider distributions of w 2 rather than v 2 , then the problem reduces to the one we already solved in the nonexpanding spacetime. Given such a normalized distribution ηðw 2 ; tÞ and the function aðtÞ, we can deduce the corresponding distribution of v 2 using χðv 2 ; tÞ ¼ dw 2 dv 2 ηðw 2 ; tÞ ¼ aðtÞ 2 ηðaðtÞ 2 v 2 ; tÞ:

In particular, the two distributions agree, as they must, when aðtÞ ≡ 1, but the time evolution of η under cosmic expansion is trivial, whereas χ does evolve according to the previous formula.

To recap, we work with a distribution over w 2 ¼ aðtÞ 2 v 2 that we call η and we can deduce the time evolution of η using ηðw

where w k ¼ aðtÞv k is the kick in w. Then, given ηðw 2 ; tÞ and the function aðtÞ, we can find the velocity distribution using Eq. (A8). The argument is now apparent and we briefly state it here. At late times, the distribution ηðw 2 ; tÞ behaves just like its static spacetime counterpart. That is to say, it develops a peak around w 2 k and has a width that is also comparable to w 2 k . The distribution of v 2 is then seen to have a peak that moves with time as v 2 k ∼ t 2=7 . The distribution of velocities, which can then be obtained by a simple change of variables, has a peak that moves as v k ∼ t 1=7 . The width of these distributions is always of order the peak position. We measure this width numerically for the velocity distribution and find width peak position ≈ 0.249;

which is the number we use in our model of the distribution function (2.4). An example of this evolution of the distribution peak and width is shown in Fig. 5.

APPENDIX B: FLUID APPROXIMATION

We describe perturbations of the graviton DM using the fluid approximation. This is adequate in our case since we are dealing with nonrelativistic DM throughout cosmic history. As such, we do not keep track of the full distribution function and its evolution but only discuss its first few moments, which amounts to using the fluid approximation. The relevant moments are defined in the standard way: ρðτÞ ¼ 4π a 4 ðτÞ Z q 2 dqEðq; τÞf 0 ðq; τÞ;

PðτÞ ¼ 4π 3a 4 ðτÞ Z q 2 dq q 2 Eðq; τÞ f 0 ðq; τÞ;

QðτÞ ¼ 4π 3a 4 ðτÞ Z q 2 dq q 4 Eðq; τÞ 3 f 0 ðq; τÞ; δρðk; τÞ ¼ 4π a 4 ðτÞ Z q 2 dqEðq; τÞf 0 ðq; τÞΔ 0 ðk; q; τÞ; δPðk; τÞ ¼ 4π 3a 4 ðτÞ Z q 2 dq q 2 Eðq; τÞ f 0 ðq; τÞΔ 0 ðk; q; τÞ ðρ þ PÞθðk; τÞ ¼ 4πk a 4 ðτÞ Z q 2 dqqf 0 ðq; τÞΔ 1 ðk; q; τÞ; ðρ þ PÞσðk; τÞ ¼ 8π 3a 4 ðτÞ Z q 2 dq q 2 Eðq; τÞ f 0 ðq; τÞΔ 2 ðk; q; τÞ ðρ þ PÞΘðk; τÞ ¼ 4πk a 4 ðτÞ Z q 2 dqq q Eðq; τÞ 2 f 0 ðq; τÞΔ 1 ðk; q; τÞ; ðρ þ PÞΣðk; τÞ ¼ 8π 3a 4 ðτÞ Z q 2 dq q 2 Eðq; τÞ q Eðq; τÞ 2 f 0 ðq; τÞΔ 2 ðk; q; τÞ; where q is the particle comoving momentum and Eðq; τÞ ≡ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi q 2 þ aðτÞ 2 m 2 p is its comoving energy. These definitions relate the multipoles of Δ to the variables of the fluid description. Using these we can derive the equations governing the fluid perturbations.

We can then proceed as usual by starting with the Boltzmann equation for df 0 =dt and perturbing it by taking f 0 → f 0 ð1 þ ΔÞ. Taking moments, we can then find the perturbation equations, and it is easy to check that these remain unchanged. They give As explained in the text, the distribution peak moves to higher velocity as ∼t 1=7 and the width also increases so that the ratio of the width to the peak position is about 0.249.

In the last equation, we have included a term called X that encapsulates the effect of truncation (i.e., Δ 3 ¼ 0) versus using the Ma-Bertschinger [52] approximation:

d ln f 0 d ln q :

In the case where a truncation is used, we have X ¼ 0; otherwise,

In addition, α ¼ 0 in the Newtonian gauge and α ¼ ð ḣ þ 6ηÞ=ð2k 2 Þ in the synchronous gauge. Finally the terms K i depend on whether we are using the Newtonian or synchronous gauge and they are given by K 0 ¼ ( φ ðNewtonianÞ;

where the metric perturbations ðh; η; ψ; ϕÞ are defined in the standard way (e.g., as in [52]). The system of equations in Eqs. (B1)-(B3) is not closed, and we close it in the same way that was discussed in [17]. We work mostly in the synchronous gauge. This is an incomplete gauge fixing and we are allowed to perform gauge transformations that are time independent. This can be a nuisance because once we get a solution, we have to check whether it can be removed by a gauge transformation. We usually resolve this ambiguity by demanding that the DM velocity perturbation θ ¼ 0. This can be confusing because all of our effects come from θ ≠ 0. How are we then to understand the use of synchronous gauge in this case?

First, let us review why we can set θ ¼ 0 in the case of CDM. The equation of motion for θ implies that θ ¼ 0 for CDM (i.e., with w ¼ 0 and δP ¼ 0 and when anisotropic stresses vanish). We then have that θð⃗ xÞ depends on ⃗

x only and we can use the remaining gauge freedom to set θ ¼ 0. This completes the gauge fixing when we have a CDM component.

In our case, we assume the presence of a very subdominant CDM component, with Ω c ∼ 10 -10 (this is in fact already implemented in CLASS). The gauge freedom is then removed by removing the velocity perturbation of this component. Because of its extremely low energy density, this fictitious component does not affect the evolution of the Universe in any way that has bearing on experimental constraints. It is then clear that the velocity perturbation of our dark gravitons is physical.

This leaves open one question: even in the presence of CDM, why not assume that there is a subdominant (also cold) component and then treat the CDM velocity perturbation as physical? The answer is that, as reviewed before, for CDM we have the equation θ ¼ 0 where we can afterwards set θ ¼ 0 by a gauge transformation. So whenever we get a nonzero constant-in-time θ for CDM, we can remove it by a gauge transformation. (Typically, there should be one linear combination of the two velocity perturbations that can be removed and another one that is physical.) What remains then is the velocity perturbation of the subdominant component which again has little bearing on cosmological evolution. In the case of dark gravitons, the θ perturbation is not constant in time and it cannot be removed by a time-independent gauge transformation. The assumption of the existence of an additional subdominant component with a vanishing velocity perturbation completely fixes the synchronous gauge.

Finally, we quote the expressions for the constant pressure, pseudopressure, and energy density used in Eqs. (3.6) and (3.7). These can be obtained by integrating the distribution function in Eq. (2.4) as in the expressions for P; ρ, and Q in the nonrelativistic limit. The result then gives p 0 , ρ 0 , and p 0 after factoring out the appropriate powers of aðtÞ and v today ðt=t age Þ 1=7 . These are

where we take η ¼ 0.25 as per the discussion in Appendix A. Here m is the DM mass, but it drops out of all expressions in the nonrelativistic limit of interest to us [see, for example, Eqs. (2.4), (3.6), and (3.7)].