Stellar cooling limits on KK gravitons
and dark dimensions

Edward Hardy, Anton Sokolov, Henry Stubbs 

Rudolf Peierls Centre for Theoretical Physics, University of Oxford,
Parks Road, Oxford OX1 3PU, UK

We revisit cooling bounds on light Kaluza-Klein (KK) gravitons, as arise in the dark dimension scenario, considering red giants, neutron stars, and supernovae. In addition to bremsstrahlung, we account for two novel production channels: resonant mixing with the in-medium photon and a pion-induced process in supernovae. The strongest limits arise from SN 1987A, with the emissivity from the pion process exceeding that from bremsstrahlung by a factor of a few. Given present uncertainties, we obtain a bound on the KK mass scale of $m_{\rm KK}\gtrsim 0.6\,{\rm eV}$ ($\gtrsim 500\,{\rm eV}$) for $2$ ($3$) extra dimensions. Improved understanding of the properties of pions in supernovae could strengthen these limits to roughly ${\rm eV}$ (${\rm keV}$). For $1$ extra dimension, the bounds are weaker than those from laboratory searches. We also show that constraints from KK graviton decays to Standard Model particles are less stringent than the cooling bounds if there is KK number violation at the level typically assumed in the dark dimension scenario, although these bounds could be strengthened by future observations.

1 Introduction

The idea that the Universe might contain relatively large extra dimensions with the Standard Model (SM) fields confined to a brane was originally motivated by the Electroweak hierarchy problem [Arkani-Hamed:1998jmv, Antoniadis:1998ig, Arkani-Hamed:1998sfv]. In this context, the fundamental higher-dimensional Planck scale $M_{D}$, which is the UV cutoff of the theory, is assumed to be close to the ${\,{\rm TeV}}$ scale. The radius $R$ of the extra dimensions must be such that the 4-dimensional Planck scale $M_{\rm Pl}$ takes its measured value, with

$$M_{\rm Pl}^{2}\simeq\left(2\pi R\right)^{n}M_{D}^{n+2}~,$$ (1)

where $n$ is the number of large extra dimensions (throughout, we assume that these are approximately toroidally compactified with a common radius $R$, but our results can easily be adapted to other setups, e.g. [Lykken:1999ms]). Although theoretically appealing, strong lower bounds on $M_{D}$ from collider searches, e.g. for signals of KK gravitons, have tempered the viability of this as a solution to the hierarchy problem.

The possibility of large extra dimensions has recently been revived in Ref. [Montero:2022prj] motivated by the value of the cosmological constant and string theory “swampland” conjectures [Vafa:2005ui]. In this proposal, known as the dark dimension, $R$ is assumed to be of order ${\rm eV}^{-1}$. For $n=1$ and $2$, the corresponding $M_{D}$ are around $10^{10}\,{\rm GeV}$ and $10^{5}\,{\rm GeV}$ respectively, which are allowed by high energy particle collider constraints, while larger $n$ are already excluded [ParticleDataGroup:2024cfk]. Crucial open questions remain concerning whether a full string theory construction of a dark dimension can be consistent with existing observational and cosmological constraints, for example on light scalar moduli fields that might typically accompany it. Moreover, further theoretical study of the swampland arguments that led to the proposal is essential [Branchina:2023ogv]. Nevertheless, given the scarcity of ideas to understand the value of the cosmological constant, it is worthwhile to analyse current constraints on dark dimensions and the prospects for future improvements or discovery.

Among the possible signals of large extra dimension scenarios, the experimental and observational implications of KK gravitons are especially interesting because these particles are unavoidable. For $n=1,2$, with $M_{D}$ well above collider scales, the laboratory bounds come from searches for deviations from Newtonian gravity at short distances [Lee:2020zjt] due to the lightest KK gravitons in the tower. These constrain $m_{\rm KK}\gtrsim 5\,{\rm meV}$, with a mild dependence on $n$ [Adelberger:2003zx]. Meanwhile, astrophysical constraints on KK gravitons have been studied extensively since the original proposal of large extra dimensions [Arkani-Hamed:1998sfv, Barger:1999jf, Cullen:1999hc, Hanhart:2000er, Hanhart:2001fx, Hannestad:2001jv, Hannestad:2001xi, Hannestad:2003yd] and give the strongest bounds for $n=2$ and $3$.

The most model-independent astrophysical bounds are due to “cooling”, in which KK gravitons produced in a hot stellar core escape, carrying away energy and thereby altering the evolution of a star (see e.g. [Raffelt:1996wa, Caputo:2024oqc] for comprehensive reviews). The strongest known cooling limits on large extra dimensions arise from the neutrino burst observed from the (presumed) protoneutron star in SN 1987A. These are weaker than the laboratory bounds for $n=1$, while for $n\geq 2$ they are stronger due to the larger number of energetically accessible KK states. Existing SN 1987A cooling constraints allow either one or two extra dimensions of the size suggested in the dark dimension scenario [Hannestad:2003yd].

Even stronger constraints arise in specific theories; for example, if the KK gravitons decay only to SM particles then observations of the heating of neutron stars lead to $m_{\rm KK}\gtrsim 3\,{\rm eV}$ for $n=2$ (the bounds for $n=1$ are again weaker than laboratory constraints) [Hannestad:2001jv, Hannestad:2001xi, Hannestad:2003yd, Casse:2003pj]. However, as noted in the original work [Hannestad:2001jv], and further analysed in [Mohapatra:2003ah, Anchordoqui:2025nmb], these limits do not apply if the KK gravitons dominantly decay to lighter KK gravitons due to a violation of KK number.

In this paper, we revisit the cooling constraints on KK gravitons from a wide range of astrophysical objects. In addition to the well-studied bremsstrahlung emission, we consider two new production processes. The first is resonant mixing of the in-medium photon (plasmon) with KK gravitons, which occurs when the dispersion relations of the two coincide. Similar production has previously led to strong limits on new light scalars from red giants and horizontal branch stars, and it is also relevant in supernovae [Hardy:2016kme, Hardy:2024gwy] (it may also play a role for white dwarfs, although we leave these to future work). The second is the production of KK gravitons in supernova from pions, specifically $\pi^{-}p\rightarrow nh$, where $h$ is the KK graviton. For axions, a similar process has been found to be slightly more important than bremsstrahlung in supernova [Carenza:2020cis], although there are uncertainties on the abundance and properties of pions in protoneutron stars.

Finally, we quantify the impact of the KK graviton decays to other KK gravitons on the model-dependent limits from decays to SM particles. This sets an upper bound on the KK number violation if cooling limits are to be dominant. Our work is complementary to the recent detailed study of Ref. [Anchordoqui:2025nmb], which analysed collider and cosmological constraints on dark dimension scenarios (we note that this reference also considered the impact of intra-tower KK graviton decays; our result for the decay rate differs from theirs, although we reach similar conclusions) and Ref. [Law-Smith:2023czn], which considers constraints from the decays of dark matter KK gravitons to the SM (focusing on the $n=1$ case). Other recent works studying the phenomenology of the dark dimension scenario include Refs. [Anchordoqui:2022txe, Anchordoqui:2022svl, Anchordoqui:2022tgp, Blumenhagen:2022zzw, Anchordoqui:2022ejw, Cui:2023wzo, Anchordoqui:2023oqm, Anchordoqui:2023tln, Anchordoqui:2023wkm, Law-Smith:2023czn, Anchordoqui:2023laz, Burgess:2023pnk, Dvali:2023zww, Anchordoqui:2024akj, Anchordoqui:2024dxu, Gendler:2024gdo, Basile:2024lcz, Casas:2024oak, Heckman:2024trz, Gross:2024zdh, Li:2024jko, Anchordoqui:2024xvl, Bedroya:2025fwh, Antoniadis:2025rck].

Conventions: we define $\kappa=\sqrt{32\pi G_{N}}$, where $G_{N}$ is Newton’s constant, $R$ is the radius of the extra dimension such that the mass of the lightest KK mode $m_{\rm KK}=1/R$, and our spin-2 polarization vectors are normalised such that $\epsilon^{s,\mu\nu}\epsilon^{s^{\prime}*}_{\mu\nu}=\delta^{ss^{\prime}}$ (which matches Ref. [Gleisberg:2003ue] and differs from version 4 of Ref. [Han:1998sg] by some factors of $2$ at intermediate steps).

The structure of this paper is as follows. In Section 2 we calculate the rate of resonant KK graviton production, which is relevant to red giants and supernovae. In Section 3, we analyse the limits from the tip of the red giant branch, and in Section 4 we obtain the bounds from old neutron stars. In Section 5 we study the bounds from SN 1987A, and in Section 6 we analyse the impact of decays within the KK graviton tower. We conclude and discuss possible future improvements in Section 7. Additional technical details are given in Appendices.

2 Thermal production and KK graviton mixing

The production rate of a weakly coupled boson $x$ in a thermal medium can be shown to be [Laine:2016hma, Hardy:2024gwy]

$$\frac{dN}{dVdt}=-g\int\frac{d^{3}k}{(2\pi)^{3}}\frac{f_{B}(\omega)}{\omega}\operatorname{Im}\left[\Pi_{xx}-\frac{(\Pi_{\gamma x})^{2}}{\Pi_{\gamma\gamma}-m_{x}^{2}}\right]+\mathcal{O}(\lambda^{4}),$$ (2)

where $\Pi_{\gamma\gamma}$ ($\Pi_{xx}$) denotes the photon’s (weakly coupled boson’s) thermal self-energy [Bellac:2011kqa, Laine:2016hma], and $\Pi_{\gamma x}$ is the 1PI mixing between them. All Lorentz indices have been contracted with appropriate polarization vectors and $g$ is the number of states for each polarization. The (assumed small) coupling between the boson and the SM is denoted $\lambda$, and $\omega=\sqrt{k^{2}+m_{x}^{2}}$.

If the mass of the new boson, $m_{x}$, exceeds the plasma frequency, $\omega_{p}$, which is the typical scale of $\Pi_{\gamma\gamma}$, the second term in the square bracket in Eq. (2) can be neglected. The first term contains all of the typical contributions; depending on the interactions of $x$, these include electron annihilation, Compton-like emission, and bremsstrahlung. This can be evaluated in various ways, for example a full finite-temperature field theory calculation using cutting rules [Kobes:1985kc, Bedaque:1996af, Gelis:1997zv], or by other systematic methods [Brahma:2025wos]. However, we will employ the simple and commonly used approach of thermally averaging the vacuum interaction rates for these processes (although this may miss interesting thermal effects such as the effective nucleon mass in supernovae).

For light beyond Standard Model (BSM) particles, the second term in Eq. (2) may be significant if the BSM particle has quantum numbers such that $\Pi_{\gamma x}$ can be non-zero. For new particles with interactions such that $\Pi_{xx}$ is related to $\Pi_{\gamma x}$ (e.g. dark photons for which $\Pi_{xx}=\lambda\Pi_{x\gamma}=\lambda^{2}\Pi_{\gamma\gamma}$) cancellations between the two terms can be important. However, this does not typically occur for more general new particle interactions (including for KK gravitons) and each term can be considered separately. Then, expanding the self-energies into real and imaginary parts, and assuming that $\operatorname{Re}[\Pi_{\gamma x}]\gg\operatorname{Im}[\Pi_{\gamma x}]$, the second term becomes

$$\frac{dN}{dVdt}=g\int\frac{d^{3}k}{(2\pi)^{3}}\frac{f_{B}(\omega)}{\omega}\operatorname{Re}[\Pi_{\gamma x}]^{2}\frac{-\operatorname{Im}[\Pi_{\gamma\gamma}]}{(\operatorname{Re}[\Pi_{\gamma\gamma}]-m_{x}^{2})^{2}-\operatorname{Im}[\Pi_{\gamma\gamma}]^{2}}.$$ (3)

Provided that the lifetime of the photon in the medium is long, such that $\operatorname{Im}[\Pi_{\gamma\gamma}]$ is small, this contribution has a sharp resonance when $m_{x}^{2}=\operatorname{Re}[\Pi_{\gamma\gamma}]$, which might dominate the total production rate.

In the narrow resonance regime, we may make a delta function approximation to give

$$\frac{dN}{dVdt}=\frac{g}{2\pi}\int d\omega\frac{k}{e^{\omega/T}-1}\operatorname{Re}[\Pi_{\gamma x}]^{2}\delta(\operatorname{Re}[\Pi_{\gamma\gamma}]-m_{x}^{2}),$$ (4)

which can be evaluated to give the energy production rate

$$\frac{dQ}{dVdt}=\frac{g}{2\pi}\frac{k_{*}\omega_{*}}{e^{\omega_{*}/T}-1}\operatorname{Re}[\Pi_{\gamma x}(\omega_{*})]^{2}\left|\frac{d\operatorname{Re}[\Pi_{\gamma\gamma}(\omega)]}{d\omega}\right|^{-1}_{\omega=\omega_{*}},$$ (5)

where $\omega_{*}$ and $k_{*}$ are the particle’s energy and momentum evaluated on-resonance. Using the standard expressions for $\operatorname{Re}[\Pi_{\gamma\gamma}]$ [Braaten:1993jw], we have

$$\left.\frac{d\operatorname{Re}[\Pi_{\gamma\gamma}^{L}(\omega)]}{d\omega}\right|_{\omega=\omega_{*}} =-\omega_{*}\frac{m_{x}^{2}}{k_{*}^{2}}\left(2+\frac{m_{x}^{2}}{\omega_{*}^{2}}\frac{(\omega_{*}^{2}-k_{*}^{2}v_{0}^{2})-3\omega_{p}^{2}}{\omega_{*}^{2}-k_{*}^{2}v_{0}^{2}}\right),$$ (6)

$$\left.\frac{d\operatorname{Re}[\Pi_{\gamma\gamma}^{T}(\omega)]}{d\omega}\right|_{\omega=\omega_{*}} =-\omega_{*}\frac{m_{x}^{2}}{k_{*}^{2}}\left(\frac{m_{x}^{2}}{\omega_{*}^{2}}-\frac{3\omega_{p}^{2}-2m_{x}^{2}}{\omega_{*}^{2}-k_{*}^{2}v_{0}^{2}}\right),$$ (7)

where $v_{0}$ is the typical velocity of electrons that contribute to the photon dispersion, equal to the Fermi velocity for degenerate plasmas.

Inside a star, $\omega_{p}$ typically decreases with increasing radial distance $r$. For longitudinal photons, a resonance always occurs if $m_{x}<\omega_{p}$, corresponding to all $r$ up to some maximum radius where $\omega_{p}(r)=m_{x}$. For transverse photons there is only a resonance if $m_{x}$ is in some narrow window above $\omega_{p}$ meaning that, if $m_{x}$ lies roughly between the smallest and largest values of $\omega_{p}$ in the star, there is a spatially thin shell where the resonance can occur. The width of this shell depends on the typical velocity of the electrons in the medium, and in the non-relativistic limit it can be treated as a delta function. This is illustrated in Figures 2 and 2, which show the dispersion relations crossing and the relevant spatial regions. The situation is more complicated in the case of KK gravitons, because there is an almost continuous spectrum of massive spin-2 bosons with both longitudinal and transverse polarization states. Both can be produced resonantly, and the range of masses means that resonant production can occur anywhere within the star.

Figure 1: The energy-momentum relation of a KK graviton of mass $m_{h}=5{\,{\rm MeV}}$ (dashed black). $m_{h}$ determines the plasma frequencies $\omega_{p}$, and hence the locations in a star, where resonant production occurs. For the longitudinal mode, there is a resonance whenever $\omega_{p}\geq m_{h}$. For transverse modes, there is a resonance when $\omega_{p}\in\left[\sqrt{2/3}m_{h},m_{h}\right]$. The regions spanned by all such longitudinal (transverse) photon dispersions in a typical protoneutron star are shown in shaded red (blue) (specifically the $20M_{\odot}$ progenitor model discussed in Section 5, which has maximum $\omega_{p}\simeq 14\,{\rm MeV}$). Representative dispersion curves are plotted in bold red (blue), and the resonance points where they intersect the KK graviton curve are marked with solid dots.

Figure 2: The regions within a protoneutron star core where the different resonances can occur for a KK graviton of mass $5{\,{\rm MeV}}$, using the same model of the star as in Figure 2. A longitudinal resonance occurs whenever $m_{x}\leq\omega_{p}(r)$, and therefore throughout the region $0<r<r_{1}$ where $\omega_{p}(r_{1})=m_{x}$ (shaded red). In the limit that the electrons are relativistic, a transverse resonance only occurs if $\omega_{p}(r)\leq m_{x}\leq\sqrt{3/2}\omega_{p}(r)$. This therefore takes place in some shell $r_{1}<r<r_{2}$ where $\omega_{p}(r_{2})=\sqrt{2/3}m_{x}$ (shaded blue). Here $r_{1}\approx 16{\,{\rm km}}$ and $r_{2}\approx 18{\,{\rm km}}$.

self_energies

(110,100) i o photoni,v1 dbl_wigglyv2,o fermion,right,tension=0.5v1,v2 fermion,right,tension=0.5v2,v1	(110,100) i o t photon,tension=100i,v1 dbl_wiggly,tension=100v1,o photonv1,v2 fermion,right,tension=0.2v2,t fermion,right,tension=0.2t,v2	(110,100) i o t1,t2,A,t3,t4 photon,tension=100i,v1 dbl_wiggly,tension=100v1,o fermionv1,t3 fermiont3,t2 fermiont2,v1
(110,60) i2,i1 o2,o1 fermioni1,v1 fermionv1,v2 fermionv2,o1 photoni2,v1 dbl_wigglyv2,o2	(110,60) i1,i2 o1,o2 photoni1,v1 dbl_wigglyv1,o1 photonv1,v2 fermioni2,v2 fermionv2,o2	(110,60) i1,i2 o1,o2 photoni1,v1 dbl_wigglyv1,o1 fermioni2,v1 fermionv1,o2
s and u channel scattering (only the s channel is shown)	t channel scattering	Four-point contact diagram

Table 1: The Feynman diagrams contributing to the photon–KK graviton mixing $\Pi_{\gamma h}$, and the cut diagrams showing the correspondence with various scattering channels. Single wavy lines are photons, and double wavy lines are gravitons. The fermion loops contain contributions from all species present in the plasma.

To determine the resonant production rate we need to evaluate $\operatorname{Re}[\Pi_{\gamma x}(\omega_{*})]$. In the case of KK gravitons, $h$, the diagrams contributing to this are shown in Table 1, together with the corresponding forward scattering processes (for the connection see, e.g., Ref. [Braaten:1993jw]).

We calculate the photon–KK graviton plasma mixing in various regimes in Appendix A. Due to charge-conjugation (C) symmetry violation being required for the mixing, this vanishes if the plasma has equal numbers of particles and antiparticles, as in the case of scalar-photon mixing [Hardy:2016kme]. However, unlike in previously studied cases, the non-relativistic mixing is suppressed by a factor of $T/m_{f}$, where $m_{f}$ is the mass of the fermion in the medium, so it vanishes in the $T\rightarrow 0$ limit. The physical justification for this can be seen from the classical equations of motion, in the limit where both the electrons and ions of the medium begin at rest. An incoming electromagnetic wave displaces the electrons and ions in opposite directions. Usually, the motion of the ions is neglected [Hardy:2016kme, Dubovsky:2015cca] because their displacement is suppressed by a factor of $m_{e}/m_{i}$. However in our case, the coupling to gravity is larger for the ions by a factor $m_{i}/m_{e}$. The two contributions to the source term for the gravitons then cancel out, leading to no mixing. The mixing is also suppressed in the limit where all particles are ultra-relativistic; this time by a factor $m_{f}^{2}/T^{2}$, and vanishes in the limit of massless fermions.

Resonant production may therefore play a role in astrophysical plasmas with a mixture of non-relativistic and relativistic particle species, provided $\omega_{p}$ is not much larger than the temperature, otherwise the Boltzmann suppression in Eq. (5) is substantial. Such systems include red giants, white dwarfs, and the cores of supernovae. For older neutron stars $\omega_{p}\gg T$, while in horizontal branch stars and in main sequence stars such as the Sun the medium is close to non-relativistic, so resonant production is negligible in these environments. For dark photons, resonant conversion is also active in the early-Universe plasma, which has been thought to play a role in cosmological constraints, although this has recently been reconsidered [Hook:2025pbn]. For KK gravitons, however, the near $C$-symmetry of this environment suppresses such processes, rendering them irrelevant.

In addition to resonant production, various processes including bremsstrahlung radiation and Compton-like emission can produce KK gravitons through the first term in Eq. (2). Their relative importance depends on the particular astrophysical system, so we defer an analysis to later sections.

3 Tip of the red giant branch

As a red giant (RG) uses up its hydrogen fuel, the mass of its degenerate helium core and hence its luminosity increases until the core becomes hot and dense enough to begin to fuse helium into carbon. Due to the degeneracy of the core, the associated increase in temperature leads to a runaway helium fusion via the triple-alpha process (the so-called helium flash), which liberates a vast amount of energy in a short period of time and ultimately causes the core to expand. The luminosity of the star then decreases, meaning that the moment of helium ignition sets the maximum luminosity of a red giant star, which is known as the tip of the red giant branch (TRGB). In the presence of any additional cooling the core would grow heavier and more luminous before the helium flash can occur, making the resulting TRGB brighter. Since the triple-alpha process has a strong dependence on the properties of the core, the TRGB brightness is a sensitive probe of new hypothetical particles that could be produced in the core and carry energy away [Raffelt:1996wa]. The main SM cooling mechanism for the core is plasmon decay to neutrinos, with an average energy loss rate just before the helium ignition of $\epsilon_{\nu}\approx 4{\,{\rm erg}}{\,{\rm g}}^{-1}\sec^{-1}$ [Straniero:2020iyi]. Recent investigations of the TRGB in Globular Clusters (GCs) indicate that the average novel energy loss can be constrained by data at the level $\epsilon_{x}\lesssim(10^{-1}\div 1)\,\epsilon_{\nu}$ [Viaux:2013lha, Capozzi:2020cbu, straniero2018axionelectroncouplingrgbtip, Straniero:2020iyi, Troitsky:2024keu, Caputo:2024oqc]. Here, we assume a conservative constraint $\epsilon_{x}<\epsilon_{\nu}$.

To estimate the rate of emission of KK gravitons from the RG core at helium ignition, we generate profiles of thermodynamic characteristics inside the core for different possible values of initial star metallicities $Z_{\text{ini}}=\left\{10^{-4},\,10^{-3},\,6\times 10^{-3}\right\}$, with the initial helium abundance $Y_{\text{ini}}=0.25$ and the age at the TRGB $\simeq\!13$ Gyr. These particular choices of parameters correspond to a recent detailed analysis of the relevant GCs in Ref. [Straniero:2020iyi]. Specifically, we evolve the one-dimensional stellar models with the MESA code [Paxton:2010ji] starting from the zero-age main sequence model with the chemical composition prescribed above and the initial mass corresponding to the above TRGB age, until the moment when helium-burning luminosity reaches $L_{\text{He}}>10^{2}L_{\odot}$. The resulting temperature, density and plasma frequency profiles are shown in Fig. 6 in Appendix B.

In the past, continuum production channels of the KK gravitons inside a RG core have been considered. The most important of these are plasmon annihilation, gravi-Compton-Primakoff scattering and gravi-bremsstrahlung. In Ref. [Barger:1999jf] the gravi-bremsstrahlung process was found to dominate over the other two. However, in this calculation the plasma was assumed to be non-degenerate, which is not a good approximation inside the highly degenerate RG core at the time of helium ignition, with such an assumption leading to an overestimation of the gravi-bremsstrahlung rate by several orders of magnitude [Hansen:2015lqa].

We find that resonant production results in an energy loss rate that is similar to the gravi-bremsstrahlung before the degeneracy correction for the latter is taken into account, and therefore we conclude that resonant emission dominates over continuum production.

Our calculations show that longitudinal resonant emission surpasses the transverse one only for lower KK masses, and that the full emission from the tower is dominated by the transverse resonance. With the resonant condition applied, the main contribution to the mixing self-energy for transverse polarizations is (see Eq. (54) in Appendix A.1):

$$\Pi_{\gamma h,f}^{T}=\frac{q_{f}e\kappa}{3\sqrt{2}\pi^{2}}p_{F}^{3}\frac{m_{h}}{v_{*}\omega_{*}}\left(1-\frac{m_{h}^{2}}{\omega_{p}^{2}}\right)\,,$$ (8)

where $f$ denotes the fermion that leads to the mixing (with $q_{f}$ its electric charge and $p_{F}$ its Fermi momentum), $m_{h}$ is the mass of a particular KK graviton mode, $\omega_{*}$ is the energy at resonance, and $v_{*}=k_{*}/\omega_{*}$ is the KK graviton velocity where $k_{*}$ is the KK graviton momentum at resonance. In a RG, the mixing is dominated by degenerate electrons, which have Fermi momentum of about $0.25\,{\rm MeV}$. Because the electrons are only mildly relativistic, the resonance occurs only for $m_{h}\simeq\omega_{p}$, as follows from the functional form of the self-energy of a transverse plasmon in degenerate plasma [Braaten:1993jw, Raffelt:1996wa]. Consequently, the factor in the brackets in the right-hand side of Eq. (8) is small $(1-m_{h}^{2}/\omega_{p}^{2})\sim 10^{-2}$, which suppresses the rate of KK graviton emission.

The resulting constraints on the size of extra dimensions are given in Table 2; evidently the assumed metallicity makes no appreciable difference. We stress that the bound on the extra energy loss that we impose, $\epsilon_{x}\lesssim 4{\,{\rm erg}}{\,{\rm g}}^{-1}\sec^{-1}$, is quite conservative. For instance, a detailed analysis for the case of axions [Straniero:2020iyi], in which axions are included in the stellar evolution, suggests $\epsilon_{a}\lesssim 0.9{\,{\rm erg}}{\,{\rm g}}^{-1}\sec^{-1}$ [Caputo:2024oqc, Carenza:2024ehj]. Moreover, a recent re-analysis [Troitsky:2024keu] using data on GCs from the Data Release 3 from the Gaia mission [Gaia:2016zol, Gaia:2023fqm] yields a further factor of 2 improvement in the constraints on the axion coupling, meaning that the energy loss bound is $\epsilon_{a}\lesssim 0.2{\,{\rm erg}}{\,{\rm g}}^{-1}\sec^{-1}$. If naively rescaled to the case of extra dimensions, this would imply an improvement to the bounds in Table 2 of a factor of $20^{1/n}$, where $n$ is the number of extra dimensions. Such a translation may be well-justified, since degenerate electrons in the RG core are good heat conductors, and so one can expect that the constraint is determined primarily by the total energy loss rather than by a detailed emission profile. Furthermore, we expect that even this constraint might be improved in the near future given the abundance of recent observational data on GCs, e.g. from the Gaia mission [Gaia:2016zol, Li_2023].

metallicity	$n=1$	$n=2$	$n=3$
$10^{-4}$	$1.0\times 10^{4}$	$2.1\times 10^{-4}$	$6.5\times 10^{-7}$
$10^{-3}$	$1.0\times 10^{4}$	$2.1\times 10^{-4}$	$6.8\times 10^{-7}$
$6\times 10^{-3}$	$1.1\times 10^{4}$	$2.2\times 10^{-4}$	$7.0\times 10^{-7}$

Table 2: Limits on the compactification radius (in meters) for $n=1,2,3$ extra dimensions for the three considered red giant profiles, labelled by the zero-age main sequence metallicity.

4 Neutron stars

Observations of the luminosities of five old neutron stars, whose ages can be inferred from kinematical data and are of order $10^{5}$ years, lead to constraints on light BSM particles when compared to predicted cooling curves. The bounds on axions are comparable to those from SN 1987A [Buschmann:2021juv], while those on new scalars can be up to four orders of magnitude stronger [Fiorillo:2025zzx]. The latter results raise hope that old neutron stars could lead to constraints on KK gravitons that are competitive with supernovae, despite the far fewer KK graviton modes that are accessible due to their much lower temperatures. In this section, we show that the resulting limits are actually somewhat weaker than from SN 1987A, although they may still provide a useful complement.

Our analysis follows that for axions of Ref. [Buschmann:2021juv] and for scalars of Ref. [Fiorillo:2025zzx] closely, so we provide only the essential inputs; further details and discussion can be found in those references. We use five neutron stars, which have observed ages $t_{\rm obs}$ and photon luminosities at infinity $L_{\rm obs}$ given in Table 3 (as discussed in detail in Ref. [Buschmann:2021juv], the luminosity measurements are more robust than surface temperature measurements that could have alternatively been used). The ages are inferred from the neutron stars’ distances from the assumed position of the supernova events.

Neutron Star	$t_{\rm obs}$/$10^{5}\,{\rm yr}$	$L_{\rm obs}$/$(10^{32}\,{\rm erg}\,{\rm s}^{-1})$
J1856	$4.2\pm 0.8$	$0.65\pm 0.15$
J1308	$5.5\pm 2.5$	$3.2\pm 0.6$
J0720	$8.5\pm 1.5$	$2.2\pm 1.2$
J1605	$4.4\pm 0.7$	$4.0\pm 1.0$
J0659	$3.5\pm 0.4$	$2.8\pm 1.4$

Table 3: The ages $t_{\rm obs}$ and photon luminosities $L_{\rm obs}$ of the five old neutron stars used in our analysis.

A cooling neutron star satisfies an energy balance equation

$$C\frac{dT_{c}}{dt}=-L_{\gamma}^{\infty}-L_{\nu}^{\infty}-L_{h}^{\infty}+H,$$ (9)

where $L_{\gamma}^{\infty}$, $L_{\nu}^{\infty}$ and $L_{h}^{\infty}$ are the photon, neutrino and graviton luminosities respectively, $T_{c}$ is the temperature of the core of the neutron star (assumed uniform throughout the core), $C$ is its heat capacity and $H$ encodes possible heating e.g. from the star’s magnetic field. In the SM, the cooling at early times is dominated by emission of neutrinos, which occurs at a rate proportional to $T_{c}^{8}$. Subsequently, surface photon emission takes over, usually when the neutron star is about $10^{5}$ years old. KK gravitons provide a new cooling channel, resulting in neutron stars that are colder than expected for their ages.

In the absence of superfluidity, the production of KK gravitons in a neutron star’s interior is expected to be dominated by bremsstrahlung (processes involving plasmons are strongly Boltzmann suppressed, except for at very early times). The emissivity of neutrino anti-neutrino pairs from neutron-proton and neutron-neutron bremsstrahlung are approximately equal in typical neutron star environments [1979ApJ...232..541F, Bottaro:2024ugp], and it is reasonable to expect the same is true for KK graviton production by bremsstrahlung. Given that (as discussed subsequently) there are uncertainties of about a factor of two on the neutrino emissivity, for our present purposes it is sufficient to make the conservative approximation of only considering KK graviton production from neutron-neutron bremsstrahlung.

We evaluate the neutron-neutron bremsstrahlung rate using the soft-radiation approximation for a degenerate neutron environment calculated in Ref. [Hanhart:2000er]. We further approximate the neutron-neutron scattering cross section, $\sigma_{0}\simeq 25\,{\rm mb}$, to be independent of energy and scattering angle. The total rate of energy emission per unit time and volume, summed over all KK graviton modes, is then given by

$$\frac{dQ}{dVdt}=\frac{\pi^{(n-17)/2}\Gamma(n+4)}{270}h_{n}Y_{n}\kappa^{2}(RT_{c})^{n}\sigma_{0}T_{c}^{4}p_{F}^{5},$$ (10)

where $p_{F}$ is the neutron Fermi momentum, and

$$h_{n} =\frac{3n^{2}+18n+19}{(n+5)(n+3)}\frac{1}{\Gamma(\frac{n+3}{2})},$$ (11)

$$Y_{n} =\zeta(n+4)+4\pi^{2}\frac{\zeta(n+2)}{(n+2)(n+3)},$$ (12)

with $\zeta$ the Riemann Zeta function. Although the degenerate and non-relativistic approximations are accurate in neutron stars, we note that Eq. (10) misses in-medium effects that, e.g., modify the neutron dispersion relation, which might be relevant (these could be interesting to analyse in the future).

We study the cooling of neutron stars using the code NSCool [2016ascl.soft09009P], which is publicly available and can be straightforwardly modified to include KK graviton emission. This code evolves a spherically symmetric model of a neutron star from early times, soon after formation, until the late times relevant for the observations of the stars that we use.

There are various nuisance parameters that complicate a comparison between theoretically predicted cooling curves and data. The most important nuisance parameters are the neutron star mass $M_{\rm NS}$, the abundance of light elements in its envelope $\Delta M$, and the equation of state in the core. An additional uncertainty is whether the neutron star’s core is superfluid, and if so what its properties are. We adopt the simple approach of assuming that there is no superfluidity and we also neglect any additional heating parametrised by $H$ in Eq. (9). These two assumptions are reasonable since they are consistent with current theoretical understanding and lead to a good fit to data in the absence of KK gravitons (as one test of this, we have checked that reversing the sign of the energy emission to KK gravitons, such that they act as a fake energy source, does not substantially improve the fit). If superfluidity was included there would be additional KK graviton production channels involving breaking and formation of Cooper pairs; these however could not lead to constraints unless it was conclusively established that superfluidity was present.

We obtain limits by evaluating the neutron star cooling curves for different values of $m_{\rm KK}$ for 6 neutron star masses linearly spaced between $M_{\odot}$ and $2M_{\odot}$, and 8 light element abundances logarithmically spaced between $10^{-20}M_{\odot}$ and $10^{-6}M_{\odot}$. We then construct a likelihood function for the measured data for each of the observed neutron stars, indexed $i$ below, where we assume the errors on the neutron star data are Gaussian,

$$\mathcal{L}^{i}(L_{\rm obs}^{i},t_{\rm obs}^{i},\sigma_{L}^{i},\sigma_{t}^{i}|m_{\rm KK},\bm{\theta})=\mathcal{N}\left(L(m_{\rm KK},\bm{\theta})-L_{\rm obs}^{i},\sigma_{L}^{i}\right)\,\mathcal{N}\left(t-t_{\rm obs}^{i},\sigma_{t}^{i}\right)~,$$ (13)

where $\mathcal{N}(x,\sigma)$ denotes the normal distribution with mean $x$ and standard deviation $\sigma$, and $\bm{\theta}=(t,M_{\rm NS},\Delta M)$ are the nuisance parameters. For this purpose, the age of the neutron star $t^{i}$ at which the predicted luminosity is evaluated is a nuisance parameter. Given a value of $m_{\rm KK}$ and set of nuisance parameters for each neutron star (which we collectively denote $\Phi\equiv({\bm{\theta}}^{i})$) the total likelihood $\mathcal{L}_{\rm tot}(m_{\rm KK},\Phi)$ is the product of the likelihoods for each neutron star.

The best fit value of $m_{\rm KK}$ (or equivalently $R$), which we call $\hat{m}_{\rm KK}$ is obtained by maximising the total likelihood, and we denote the associated values of the nuisance parameters $\hat{\Phi}$. We obtain a lower bound on $m_{\rm KK}$ from the test statistic

$$q(m_{\rm KK})=-2\ln\left(\frac{\mathcal{L}_{\rm tot}(m_{\rm KK},\Phi_{\rm max}(m_{\rm KK}))}{\mathcal{L}_{\rm tot}(\hat{m}_{\rm KK},\hat{\Phi})}\right),$$ (14)

where $\Phi_{\rm max}(m_{\rm KK})$ are the nuisance parameters that maximise the likelihood for a given $m_{\rm KK}$. By Wilks’ theorem, the $95\%$ confidence lower bound on $m_{\rm KK}$ is given by $q(m_{\rm KK})=2.7$ (the non-applicability of Wilks’ theorem might lead to an order-one error in our limits, which is not a concern for our present purposes [Buschmann:2021juv]).

We carry out the procedure above for three different choices of equation of state, given the uncertainty on this: the models APR [Akmal:1998cf], BSk22, and BSk26 [Pearson:2018tkr]. We find that the BSk22 equation of state leads to slightly weaker limits on $m_{\rm KK}$ than the others (by a factor of about $1.5$ for $n=2$ extra dimensions), and it fits the data slightly better in the absence of KK gravitons, so we use this for the limits that we quote.

$n=1$	$n=2$	$n=3$
$1.5\times 10^{3}$	$7.5\times 10^{-5}$	$2.7\times 10^{-7}$

Table 4: Limits on the compactification radius (in meters) for $n=1,2,3$ extra dimensions from the five old neutron stars that we analyse.

The $95\%$ confidence level bounds obtained this way are given in Table 4. As an indication of the impact of KK gravitons on the evolution of a neutron star, in Figure 3 we plot the best fit cooling curves for the neutron star J1605 in the absence of KK gravitons and for $n=2$ extra dimensions with $R=2\times 10^{-4}\,{\rm m}$ (i.e. $m_{\rm KK}=10^{-3}\,{\rm eV}$, slightly excluded by our constraint). In both these cases, the best fit curve has a neutron star mass of $M_{\odot}$ and light element abundance of $10^{-12}M_{\odot}$. We also show the energy emission to photons, neutrinos and KK gravitons as a function of time. KK graviton emission dominates between approximately $10^{4}$ and $10^{5}$ years, around when the dominant SM emission changes from neutrinos to surface photons.

Figure 3: Left: The observed photon luminosity and age of the neutron star $J1605$ with uncertainties (red), and the best fit cooling curves assuming only Standard Model energy loss (blue, “SM”) and including also energy loss to KK gravitons with $2$ extra dimensions of radius $2\times 10^{-4}\,{\rm m}$ (dashed purple). Combining the likelihoods over all $5$ neutron stars, the latter theory is excluded. Right: The energy loss to photons ($L_{\gamma}^{\infty}$), neutrinos ($L_{\nu}^{\infty}$), and KK gravitons ($L_{h}^{\infty}$) in the same large extra dimensional theory as plotted in the left panel.

We note that it might be interesting to calculate the production rate in the one-pion-exchange approximation [Cullen:1999hc] with in-medium corrections. This would be fully consistent with the way that neutrino production from bremsstrahlung is treated in NSCool, which builds on the expressions obtained in Ref. [1979ApJ...232..541F] (Refs. [Buschmann:2021juv] and [Fiorillo:2025zzx] carry out analyses along these lines for axions and scalars). As mentioned, it would also be useful to include KK graviton production by neutron-proton bremsstrahlung. Given the uncertainties on the equation of state, the fact that the neutrino emission rate likely differs from that implemented in NSCool by factors of about two [Bottaro:2024ugp], and that the limits we obtain are somewhat weaker than those from SN 1987A, we do not regard these improvements as essential for our present work.

Constraints on the couplings of axions have also been obtained from observations of young neutron stars, such as Cassiopeia A (Cas A) [Leinson:2014ioa, Sedrakian:2015krq, Hamaguchi:2018oqw, Bhat:2019tnz, Leinson:2021ety], which turn out to be similar to the bounds from old neutron stars. Since the rate of energy emission to axions scales as $T_{c}^{6}$, from Eq. (10) we expect that for $n=1$ and $2$ the limits on KK gravitons from young neutron stars will be comparable to those we have obtained from old stars. For larger $n$ the limits from young neutron stars might be somewhat stronger. However, since the bounds from supernova are much more important for large $n$, we do not consider young neutron stars here (there are presently also uncertainties on the Cas A observational data due to instrumentation issues with the Chandra telescope [Posselt:2018xaf, Buschmann:2021juv]).

5 SN 1987A

During a core-collapse supernova, neutrinos are produced abundantly over a cooling period lasting approximately ten seconds. The high temperatures (up to roughly $50\,{\rm MeV}$) and densities (of order $10^{14}\,{\rm g}\,{\rm cm}^{-3}$) of the protoneutron star during this time make it an ideal environment for the production of low-mass BSM particles. Such particles would lead to additional cooling, shortening the duration of the neutrino cooling era. Despite limited observational data and relatively imprecise theoretical understanding of the dynamics of the explosion, the efficiency of BSM particle production means that supernova constraints are often competitive with other astrophysical bounds, and moreover provide the only cooling limits at higher masses. Even though the relevant values of $m_{\rm KK}$ are small, the high supernova temperatures provide a significant advantage in constraining the large-extra-dimension scenario, because many KK modes are kinematically accessible.

The key observational data comes from SN 1987A, which, to date, is the only core-collapse supernova for which a neutrino burst was detected and can hence be compared to theoretical predictions. Raffelt proposed an analytic constraint that the new energy loss should not exceed $10^{19}{\,{\rm erg}}{\,{\rm g}}^{-1}\sec^{-1}$ [Raffelt:1987yt, Raffelt:1990yz]. This criterion has since been widely applied, including to various theories of new scalars [Krnjaic:2015mbs, Balaji:2022noj, Hardy:2024gwy], axions (see [Caputo:2024oqc] and references therein), dark photons (e.g. in [Chang:2016ntp]), as well as KK gravitons [Cullen:1999hc, Hanhart:2000er, Hanhart:2001fx, Hannestad:2003yd].

However, significant uncertainties remain. The explosion mechanism of SN 1987A has not yet been confirmed; although observations are consistent with a remnant neutron star (as expected following core-collapse), other possibilities are not robustly excluded, and in such cases the cooling bound may not apply [Bar:2019ifz]. Even if SN 1987A was a standard core-collapse event, the physics of the explosion is complex and numerical simulations still do not fully capture all processes. As a result, the energy-loss criterion that we use is only reliable to within factors of order-one (indeed slightly different versions of it are used in the literature).

We calculate KK graviton production using supernova profiles from three representative spherically symmetric numerical simulations [Rampp:2002bq, Woosley:2002zz, Buras:2005rp] available from the Garching Core-Collapse Supernova Archive [GarchingArchive], labelled by their progenitor mass (further discussion of these can be found in Ref. [Hardy:2024gwy]). We take the profiles immediately prior to when the shock revival is expected to occur and consider energy loss out to the neutrino sphere, at a radius of about $30\,{\rm km}$, as identified in the simulations, imposing the Raffelt criteria in the form above over this volume. The typical core temperatures and densities in the profiles used are consistent with those found in modern three-dimensional simulations during the first few seconds after bounce [Burrows:2024pur].

We consider three production mechanisms: the well-studied bremsstrahlung channel, as well as the new resonant production and pion-induced processes. The resulting upper bounds on the compactification radius, $R$, are summarised in Table 5 (the assumptions underlying the two quoted $\pi^{-}$ process limits are discussed below). In this table, we also show the previous constraints of Refs. [Hannestad:2003yd] and [Hanhart:2000er], which consider bremsstrahlung production using a constant temperature and density supernova model. We see that the limits from the pion process are slightly stronger than those from bremsstrahlung, while those from resonant production are much weaker.

		$n=1$	$n=2$	$n=3$
	$17.8M_{\odot}$      Brem.	$3.0\times 10^{2}$	$5.7\times 10^{-7}$	$6.6\times 10^{-10}$
	Res.	$9.6\times 10^{4}$	$2.3\times 10^{-5}$	$1.6\times 10^{-8}$
	$\pi^{-}$ (a)	$3.6\times 10^{2}$	$3.8\times 10^{-7}$	$4.1\times 10^{-10}$
	$\pi^{-}$ (b)	$2.3\times 10^{2}$	$2.7\times 10^{-7}$	$2.6\times 10^{-10}$
	$20M_{\odot}$      Brem.	$5.3\times 10^{2}$	$8.0\times 10^{-7}$	$8.5\times 10^{-10}$
	Res.	$1.3\times 10^{5}$	$2.7\times 10^{-5}$	$1.8\times 10^{-8}$
	$\pi^{-}$ (a)	$7.8\times 10^{2}$	$5.8\times 10^{-7}$	$5.5\times 10^{-10}$
	$\pi^{-}$ (b)	$5.4\times 10^{2}$	$4.6\times 10^{-7}$	$4.2\times 10^{-10}$
	$25M_{\odot}$      Brem.	$8.3\times 10^{1}$	$2.7\times 10^{-7}$	$3.7\times 10^{-10}$
	Res.	$4.6\times 10^{4}$	$1.6\times 10^{-5}$	$1.2\times 10^{-8}$
	$\pi^{-}$ (a)	$5.6\times 10^{1}$	$1.4\times 10^{-7}$	$2.0\times 10^{-10}$
	$\pi^{-}$ (b)	$3.8\times 10^{1}$	$1.2\times 10^{-7}$	$1.8\times 10^{-10}$
	Hannestad & Raffelt [Hannestad:2003yd]	$4.9\times 10^{2}$	$1.0\times 10^{-6}$	$1.1\times 10^{-9}$
	Hanhart et al. [Hanhart:2001fx]		$7.1\times 10^{-7}$	$8.5\times 10^{-10}$

Table 5: The limits on the compactification radius (in meters) for $n=1,2,3$ extra dimensions for three example supernova models, labelled by their progenitor mass $17.8M_{\odot}$, $20M_{\odot}$ and $25M_{\odot}$ (corresponding to protoneutron star masses of $1.59M_{\odot}$, $1.50M_{\odot}$ and $1.86M_{\odot}$ respectively). “Brem.” indicates the limit using only nucleon-nucleon bremsstrahlung emission, and “Res.” indicates the limit using only resonant emission. $\pi^{-}$ (a) and (b) refer to limits considering only the pion process $\pi^{-}p\rightarrow nh$ with the two sets of approximations described in the main text.

Taking the $17.8M_{\odot}$ progenitor, which yields the weakest bounds, and combining emission from bremsstrahlung and pions with the more conservative assumptions (a), leads to limits of $m_{\rm KK}\gtrsim 0.6\,{\rm eV}$ for $n=2$ and $m_{\rm KK}\gtrsim 500\,{\rm eV}$ for $n=3$. These are somewhat stronger than those of Refs.[Hanhart:2001fx] and [Hannestad:2003yd], but, given the substantial uncertainties, the differences are not dramatic. The stronger limits arising from the $25M_{\odot}$ progenitor suggest that improved theoretical understanding of the explosion, and better understanding of the properties of pions in supernovae, might allow stronger limits to be placed. Should we be lucky enough to witness another supernova in our Galaxy or a nearby dwarf galaxy, the resulting high-precision neutrino observations would certainly also be valuable.

5.1 Bremsstrahlung

We evaluate the bremsstrahlung production rate following Refs. [Hanhart:2000ae, Hanhart:2000er] (see also [Hannestad:2003yd]) using the soft-radiation approximation, which is an expansion in the parameter $\chi\equiv\omega m_{N}/p^{2}$, where $\omega$ is the energy of the emitted KK graviton and $p^{2}/m_{N}$ is the energy of the nucleon-nucleon system (with $p$ the centre of mass momentum and $m_{N}$ the nucleon mass). This approximation amounts to considering only the emission from the external nucleon legs, which dominates at small $\chi$ (with rate scaling as $\chi^{-2}$). In general there will be many contributions at order $\chi^{0}$, which we neglect. We also note that this approach misses multi-body and in-medium correlation effects.

We further assume that the nucleon-nucleon scattering cross section is energy independent and has no dependence on the scattering angle, and that nucleons in the protoneutron star can be approximated as non-degenerate and consist solely of neutrons. We expect these approximations to lead to at most corrections of order unity to the rate. The rate of energy emission to KK gravitons by bremsstrahlung is then given by [Hanhart:2000er]

$$\frac{dQ}{dVdt}=\frac{\kappa^{2}}{18}\sqrt{\frac{T^{7}}{m_{N}}}(RT)^{n}\pi^{\frac{n-3}{2}}h_{n}n_{n}^{2}\sigma_{0}^{nn}\frac{\Gamma(n+1)\Gamma(n+5)}{\Gamma(n+\frac{5}{2})}\frac{n^{2}+9n+23}{(2n+7)(2n+5)},$$ (15)

where $h_{n}$ is defined in Eq. (11), $n_{n}$ is the neutron number density, $\sigma_{0}^{nn}\approx 25\text{mb}$, and we approximate $m_{N}$ by its vacuum value.

As shown in [Hanhart:2000er], at a temperature of $30\,{\rm MeV}$ the expansion parameter $\chi\simeq 3/4$, so in the hottest regions of protoneutron stars (with temperatures around $40\,{\rm MeV}$) the soft approximation is not fully reliable. This is especially the case because the most massive KK modes are the most numerous for $n>1$. Nevertheless, we regard it as a more reliable and conservative approach than a simple one-pion-exchange calculation, since detailed studies indicate that the latter overestimates the emission rate of axions by up to about an order of magnitude [Chang:2018rso, Carenza:2019pxu] (because, e.g., the corrections from multi-pion exchange are substantial). For the case of KK graviton emission, the one-pion-exchange result also leads to an energy emission rate almost ten times larger than the soft-approximation result [Hanhart:2000er]. A future detailed analysis of the accuracy of the soft approximation would be useful, but is beyond the scope of our present work.

5.2 Resonant production

The rate of KK graviton production by resonant mixing is calculated as described in Sections 2 and Appendix A. For the profiles that we consider, the plasma frequency takes values up to roughly $16{\,{\rm MeV}}$, which hence sets the maximum KK graviton mass that can be produced by this process. We find that, as for the case of red giants, the longitudinal resonance only dominates for lighter KK modes, and, when summed over the full tower, the transverse resonance is responsible for almost all of the energy production. This effect is more pronounced at larger $n$ since the heavier KK modes have a larger multiplicity.

As can be seen from Table 5, the total rate of energy emission by resonant production is small compared to bremsstrahlung, in contrast to what is found for, e.g., dark photons. We attribute this to the combination of several factors. First, the number of KK modes that are accessible for bremsstrahlung goes as $(RT)^{n}$, whereas in the resonant process the number scales as $(R\omega_{p})^{n}$, which suppresses the total emission in the regions where $\omega_{p}\ll T$. In a supernova core, there are regions where $\omega_{p}\sim T$ due to the large electron chemical potential; however, due to the shape of the photon dispersion, the resonant frequency $\omega_{*}$ is always greater than $\omega_{p}$ and hence $T$ by a potentially large factor, unless the graviton mass $m_{h}\sim\omega_{p}$. This incurs a large Boltzmann suppression factor $e^{-\omega_{*}/T}$. Additionally, even focussing on a specific KK mode with $m_{h}\sim\omega_{p}$ the expression for the transverse mixing contains a factor of $(1-m_{h}^{2}/\omega_{p}^{2})$ which suppresses the production, and similarly the longitudinal production is suppressed (see Eqs. (55) and (56) in Appendix A.1). Numerically we find that removing the Boltzmann suppression factor as well as the factor that suppresses the mixing increases the production rate from transverse resonance by a factor of roughly $300$ for $n=2$, explaining most of the observed difference between the resonant and bremsstrahlung limits.

We also note that an interesting related process is the annihilation of two plasmons into a KK graviton. Previous calculations rely on a phenomenological photon mass term [Das:2008ss], which violates a Ward identity that we would expect in a consistent thermal field theory calculation, or neglect the plasma mass of the photon [Barger:1999jf]. The resulting SN 1987A bounds are somewhat weaker than those arising from bremsstrahlung, but are free from nuclear physics uncertainties. It would be beneficial to revisit this by evaluating the $\gamma\gamma h$ vertex consistently within thermal field theory.

5.3 Production from pions

In a cooling protoneutron star, the number density of negatively charged pions $n_{\pi^{-}}$ is typically a few percent of that of nucleons [Fore:2019wib]. Such a pion number density is present, even though the temperature is well below the pion mass, because there is a substantial isospin asymmetry. This leads to a chemical potential of negative pions $\mu_{\pi^{-}}=\mu_{n}-\mu_{p}$ (enforced by $n\leftrightarrow p\pi^{-}$ equilibrium) that is typically well above $50\,{\rm MeV}$ and may even approach $m_{\pi}$. In contrast, the number density of positive and neutral pions is negligible. It has recently been shown that the number density of negative pions can be further enhanced by attractive interactions with nucleons in the medium, which reduce the effective pion energy [Fore:2019wib].

The negative $\pi^{-}$ abundance in supernovae can lead to axion emission that is comparable to that from bremsstrahlung through $p\pi^{-}\rightarrow na$; see [Turner:1991ax, Raffelt:1993ix, Mayle:1987as] for early work of this process and [Carenza:2020cis] for a recent study demonstrating its significance (similar production of dark photons has been considered in Ref. [Shin:2022ulh]). Here we analyse the analogous channel $p\pi^{-}\rightarrow nh$. It is plausible that such a process might be relatively more important for KK gravitons than axions. This is because the emission spectrum is expected to be harder than bremsstrahlung due to the pion energy, which opens up additional channels in the case of KK gravitons (a recent update to the classic paper [Hannestad:2003yd] highlighted the pion-induced process as potentially relevant for this reason).

At leading order in chiral perturbation theory, the vacuum matrix element for $p\pi^{-}\rightarrow nh$ is given by the sum of four terms $M_{\alpha\beta}=M^{s}_{\alpha\beta}+M^{u}_{\alpha\beta}+M^{t}_{\alpha\beta}+M^{c}_{\alpha\beta}$ where $M^{c}_{\alpha\beta}$ is the 4-point $np\pi^{-}h$ interaction, corresponding to the Feynman diagrams in Figure 4. The Feynman rules required to evaluate these are given in Ref. [Han:1998sg], except for the 4-point vertex, which can be derived from the chiral perturbation theory Lagrangian (see Eq. (60) in Appendix C).

Pion_diagrams

(90,60) i2,i1 o2,o1

label=$\pi^{-}$i1 label=$p$i2 label=$h$o1 label=$n$o2

fermion,label=$p_{2}$i2,v1 fermionv1,v2 fermion,lab.side=right,label=$p_{4}$v2,o2 dashes,label=$p_{1}$i1,v1 dbl_wiggly,lab.side=left,label=$p_{3}$v2,o1

(a)

(90,60) i2,i1 o2,o1

label=$\pi^{-}$i1 label=$p$i2 label=$h$o1 label=$n$o2

fermion,label=$p_{2}$i2,v1 fermionv1,v2 fermion,lab.side=right,label=$p_{4}$v2,o2

phantomi1,v1 phantomv2,o1

dashes,label=$p_{1}$,lab.side=lefti1,u1 dashesu1,v2 dbl_wigglyv1,u2 dbl_wiggly,lab.side=left,label=$p_{3}$u2,o1

(b)

(90,60) i2,i1 o2,o1

label=$\pi^{-}$i1 label=$p$i2 label=$h$o1 label=$n$o2

fermion,label=$p_{2}$,lab.side=lefti2,v1 fermion,label=$p_{4}$,lab.side=leftv1,o2 dashesv1,v2 dashes,label=$p_{1}$i1,v2 dbl_wiggly,label=$p_{3}$,lab.side=rightv2,o1

(c)

(90,60) i2,i1 o2,o1

label=$\pi^{-}$i1 label=$p$i2 label=$h$o1 label=$n$o2

fermion,label=$p_{2}$,label.side=lefti2,v1 fermion,label=$p_{4}$,label.side=leftv1,o2 dashes,label=$p_{1}$,label.side=lefti1,v1 dbl_wiggly,label=$p_{3}$,label.side=leftv1,o1

(d)

Figure 4: The four Feynman diagrams contributing to the pion process $p\pi^{-}\rightarrow nh$.

The expression for the matrix element is given in Appendix C. The total matrix element obeys a Ward identity $p_{3}^{\alpha}M_{\alpha\beta}=p_{3}^{\beta}M_{\alpha\beta}=0$, which guarantees that terms which diverge as $m_{h}\rightarrow 0$ cancel (this requires that the internal nucleons are not approximated as on-shell). The fermion spin sum can then be performed, as well as sums over graviton polarizations as described in Ref. [Han:1998sg]. The resulting full expression is long, and so we take the large $m_{N}$ limit to give

$$\sum|M|^{2}=\frac{\kappa^{2}f_{pp}^{2}m_{N}^{2}}{6(x-1)^{2}} [4(6y-x(6x^{2}-6xy+y(y+4)+10))-4z^{2}(1-x+y)^{2}$$ (16)

$$- z(5x^{3}+x^{2}(13-8y)+x(y^{2}+12y-9)+2y^{3}-13y^{2}+12y-9)],$$

where

$$x=\frac{t}{m_{\pi}^{2}},\quad y=\frac{m_{h}^{2}}{m_{\pi}^{2}},\quad z=\frac{t-m_{\pi}^{2}}{E_{\pi}^{2}},$$ (17)

$t=(p_{4}-p_{2})^{2}$ is the usual Mandelstam variable, with the conventions for the momenta as in Figure 4, and $E_{\pi}$ is the pion energy (the unusual variables $x$, $y$ and $z$ are introduced for compactness). Numerically, this large $m_{N}$ approximation reproduces the full spin summed matrix element up to at worst $\mathcal{O}(1)$ differences.

The production rate of a single KK graviton mode per unit volume and energy $\omega$ is given by

$$\frac{d\dot{n}}{d\omega}=\int\frac{d^{3}p_{2}}{(2\pi)^{3}2m_{N}}\frac{d^{3}p_{1}}{(2\pi)^{3}2E_{\pi}}\frac{d^{3}p_{4}}{(2\pi)^{3}2m_{N}}\frac{4\pi\omega^{2}}{(2\pi)^{3}2\omega}(2\pi)^{4}\delta^{4}(p_{1}+p_{2}-p_{3}-p_{4})|M|^{2}f_{p}f_{\pi}(1-f_{n})~,$$ (18)

where $f_{i}$ are the Fermi-Dirac and Bose-Einstein distributions as appropriate, and $|M|^{2}$ is spin summed rather than averaged. Eq. (18) can be partially evaluated (analogously to the calculation in Ref. [Carenza:2020cis]), such that the total energy emission rate to a single KK mode per unit time and volume $\frac{dQ}{dVdt}=\int d\omega\,\omega\frac{d\dot{n}}{d\omega}$ is given by

$$\frac{dQ}{dVdt}=\frac{T^{11/2}}{(2\pi)^{5}\sqrt{2m_{N}}}\left(\frac{z_{p}z_{\pi}}{z_{n}+1}\right)\left(\int_{0}^{\infty}dx_{p}\frac{x_{p}}{e^{x_{p}^{2}}+z_{p}}\right)\left(\int_{x_{\rm min}}^{\infty}dx_{\pi}\frac{x_{\pi}^{2}\epsilon_{\pi}}{e^{\epsilon_{\pi}-y_{\pi}}-z_{\pi}}\int_{-1}^{1}d\cos\theta\,|M|^{2}\right),$$ (19)

where the fugacities $z_{i}=\exp((\mu_{i}-m_{i})/T)$; the dimensionless variables $x_{p}=p_{2}/\sqrt{2m_{N}T}$, $x_{\pi}=p_{1}/T$, $y_{\pi}=m_{\pi}/T$, $\epsilon_{\pi}=E_{\pi}/T$; and $x_{\rm min}$ is set by the minimum pion momentum such that the process is kinematically allowed. $\theta$ denotes the angle between the incoming pion and outgoing graviton in the centre-of-momentum frame, and the kinematics are expanded in small $m_{\pi}/m_{N}$, so $t$ is a function only of $E_{\pi}$ and the masses. A sum over KK modes leads to the total energy emission.

The most robust approach to obtain limits would be to either incorporate KK graviton production directly into a simulation that includes pions, or to evaluate Eq. (19) using a profile from a simulation that include pions, adopting the same pion and nucleon equations of state and dispersion relations as used in the simulation. The matrix element should likewise be evaluated including consistent in-medium corrections.

Unfortunately, this ideal is currently unachievable for several reasons. First, calculating in-medium corrections to the matrix element presents a major challenge. Even for axion emission via bremsstrahlung this is difficult, and KK graviton emission from pions is likely to be even more challenging. Since $\mathcal{O}(1)$ corrections to the leading-order chiral-perturbation-theory result are likely, this limits the accuracy that can realistically be reached. Additionally, even state-of-the-art simulations do not yet include pions.

Given this situation, we adopt an approach expected to provide a reasonable approximation to a full analysis. We use the proton and neutron number-density and temperature profiles from the spherically symmetric simulations discussed above, and we take the pion chemical potential to be $\mu_{\pi^{-}}=\mu_{n}-\mu_{p}$, using the neutron and proton chemical potentials determined in the simulations. For simplicity we evaluate Eq. (19) using the free-particle dispersion relations for the protons and neutrons. We fix the proton and neutron fugacities $z_{n,p}$ to reproduce their number densities as extracted from the simulations with this dispersion relation.

Given the previously discussed medium-induced enhancement of the pion number density, we present results using two approaches to the pion abundance, both using $\mu_{\pi}$ as extracted above, labelled (a) and (b) in Table 5. In the first, we use the vacuum pion dispersion relation. In the second, following Ref. [Fore:2019wib], we take the pion to have a modified dispersion relation

$$E_{\pi}(p)=\sqrt{p^{2}+m_{\pi}^{2}}+\Sigma_{\pi^{-}}(p)~,$$ (20)

where the self-energy is dominated by interactions with neutrons and is estimated from the one-loop expression

$$\Sigma_{\pi^{-}}(p)=\int\frac{d^{3}k}{(2\pi)^{3}}f_{n}(E_{n}(k))V_{n\pi^{-}}(p_{\rm cm})~,$$ (21)

with $E_{n}(k)$ the energy of a neutron of momentum $k$, and $p_{\rm cm}$ is the centre-of-mass momentum $p_{\rm cm}=\tilde{m}\left(p^{2}/m_{\pi}^{2}+k^{2}/m_{N}^{2}-2{\bf p}.{\bf k}/(m_{\pi}m_{N})\right)^{1/2}$ with $\tilde{m}$ the pion-neutron reduced mass. The pion–neutron potential is obtained from experimentally measured scattering phase shifts $\delta_{l,\nu}$ as

$$V_{n\pi^{-}}(p_{\rm cm})=-\alpha\sum_{l,\nu}(2l+1)\frac{2\pi\delta_{l,\nu}}{\tilde{m}p_{\rm cm}}~,$$ (22)

where $l$ and $\nu$ are, respectively, the angular momentum and nucleon spin projection. The factor $\alpha\simeq 0.16$ is phenomenological and is chosen such that, in the limit of small fugacities, the pion number density matches that obtained from a virial expansion. Further details can be found in Ref. [Fore:2019wib]. $\Sigma_{\pi^{-}}$ reduces the pion energy, especially at momenta around $300\,{\rm MeV}$, which enhances the pion abundance for a given $\mu_{\pi}$.

Additionally, we exclude the central regions of the supernova profiles where the pion chemical potential exceeds $120\,{\rm MeV}$ because the pions are close to condensation in these regions, and both the vacuum dispersion relation and Eq. (20) are expected to be unreliable.

The first of these two approaches is expected to be conservative, whereas the second is an uncontrolled extrapolation, although it appears to agree well with data in the regimes where comparisons are possible. In combination, they are expected to offer a reasonable estimate of the uncertainties and of the range in which the result from a full analysis would likely lie. Given the uncertainties, the limits we quote in the abstract and the start of this section are obtained from the first approach (see Appendix F of Ref. [Fiorillo:2025gnd] for a detailed useful discussion, and for Ref. [Fore:2023gwv] recent related work). As a further cross-check, we have taken the total baryon density and temperature profiles from simulations and computed the neutron, proton, and pion abundances predicted by the virial-inspired model (imposing charge neutrality and beta equilibrium), as in Ref. [Fore:2019wib]. This yields limits that are consistent with the second approach.

The results in Table 5 show that the corrected pion dispersion relation leads to a moderate strengthening of the bounds. However, the difference between the profiles is more important, primarily due to the higher temperature of the $25M_{\odot}$ progenitor. Consequently, in the future it may be useful to revisit our approach of evaluating the emission only at a single moment rather than through the first few seconds of cooling. This would require making use of profiles from three-dimensional simulations, since the spherically symmetric ones we utilise do not successfully produce explosions.

As anticipated, the energy spectrum of KK gravitons produced by the pion process is harder than that from bremsstrahlung. To illustrate this, in Figure 5 we plot the energy emission spectrum by bremsstrahlung and the pion process for a single KK mode of mass $50\,{\rm MeV}$ at a typical point in the protoneutron star. In particular, we set the temperature to $40\,{\rm MeV}$, the baryon number density to $0.1\,{\rm fm}^{-3}$, the proton fraction to $0.3$, and the pion chemical potential to $80\,{\rm MeV}$. The spectrum from the pion process has support only for KK graviton energies greater than the pion mass, because the majority of the incoming pion’s energy goes to the KK graviton. For the single KK mode shown, the emissivity from bremsstrahlung and the pion process are similar, but the total energy from the latter is greater because its harder spectrum allows more KK modes to be accessed. Relatedly, the pion process produces, on average, heavier KK gravitons than bremsstrahlung.

Figure 5: The energy emission spectrum $dQ/(d\omega dVdt)$ to a KK graviton of mass $50\,{\rm MeV}$ (with energy $\omega$) at a typical point in a protoneutron star (see main text). We show the spectra from bremsstrahlung emission (“Brem.”) and the pion process. For the latter, we plot the spectra obtained with the two sets of assumptions described in the main text (labelled “$\pi^{-}$(a)” and “$\pi^{-}$(b)”).

6 Signals from decays

In addition to the minimal cooling bounds discussed so far, constraints can arise if some of the KK gravitons produced in an astrophysical object decay to SM particles. The decay rate of a KK graviton of mass $m_{h}$ to approximately massless SM particles is typically

$$\Gamma^{\rm SM}=\lambda^{2}\frac{m_{h}^{3}}{M_{\rm Pl}^{2}}\simeq(10^{14}\,{\rm yr})^{-1}\lambda^{2}\left(\frac{m_{h}}{{\rm MeV}}\right)^{3}~,$$ (23)

where the parameter $\lambda$ encodes SM coupling constants and loop factors, so is not much smaller than one. Given Eq. (23), only the relatively heavy KK gravitons that are produced in supernovae have decay rates fast enough to yield observable signals, while the lighter modes from other types of stars decay far too slowly to be detectable.

Constraints from KK graviton decays to SM particles are strongest when no other decay channels are available. For this to be the case, decays into lighter KK gravitons must be forbidden; whether this is the case is model-dependent. In the limit of an exactly toroidal compactification, such decays are prevented by KK number conservation, the discrete remnant of momentum in the extra dimensions (there is one KK number associated with each extra dimension). In this scenario, constraints arise from various observations associated with supernovae [Hannestad:2003yd]: searches for the products of decays of KK gravitons gravitationally retained by neutron stars (such as RX J185635–3754); searches for the products of decays from the vicinity of the recent supernova remnant Cas A; searches for the products of decays from all past supernovae; and from bounds on the heating of old neutron stars by the products of gravitationally bound KK graviton decays. For $n=2$, these different signals lead to limits on $m_{\rm KK}$ that are a factor of a few to ten stronger than that from the cooling of SN 1987A [Hannestad:2003yd].

However, as discussed in [Gonzalo:2022jac, Law-Smith:2023czn], the details of the compactification may induce a small explicit breaking of KK number in the low-energy effective theory. In this case, some of the heavy KK gravitons produced in a supernova would decay to lighter KK gravitons. These lighter particles can still decay to SM particles, but at a much slower rate due to the $m_{h}^{3}$ dependence in Eq. (23), resulting in weaker constraints. It is therefore interesting to ask whether, with the magnitude of KK number breaking usually assumed in the dark dimension scenario [Gonzalo:2022jac], constraints from decays into SM particles are relevant.

We follow Ref. [Gonzalo:2022jac] in supposing that the details of the compactification allow decays to violate KK number by up to an amount that we call $\delta n$. This is assumed to be much smaller than the KK numbers of typical gravitons produced in a supernova, which is of order $10\,{\rm MeV}/m_{\rm KK}$. We label the gravitons with KK numbers $(l_{1},l_{2},\ldots),(m_{1},m_{2},\ldots),(n_{1},n_{2},\ldots)$, where the indices correspond to the extra dimensions and $l_{i},m_{i},n_{i}\in\mathbb{Z}$. In the absence of KK number violation, the three-point vertex involving KK, is non-zero only if

$$l_{i}=m_{i}+n_{i}\quad{\rm for~all~}i~.$$ (24)

(in the alternative convention of considering only positive KK numbers the corresponding condition is $l_{i}=|m_{i}\pm n_{i}|$). We therefore assume that a KK graviton $(l_{i})$ can decay to modes $(m_{i}),(n_{i})$ provided

$$l_{i}-\delta n<n_{i}+m_{i}<l_{i}+\delta n\quad{\rm for~all~}i~.$$ (25)

Again following [Gonzalo:2022jac], we suppose that the typical rate for such a decay of a KK graviton of mass $m$ is

$$\Gamma\simeq\beta^{2}\frac{m_{h}^{3}}{M_{\rm Pl}^{2}}\sqrt{\frac{m_{\rm KK}\delta n}{m_{h}}}~,$$ (26)

where $\beta$ is an unknown numerical parameter expected to be $\mathcal{O}(1)$. The factor $\sqrt{m_{\rm KK}\delta n/m_{h}}$ in Eq. (26), which is the typical velocity of the products of the decays, arises from the phase space since the decays are near threshold.

We focus on the case of two extra dimensions, because for a single extra dimension the astrophysical bounds are much weaker than laboratory constraints, even if KK gravitons decay exclusively to SM particles [Hannestad:2003yd]. We assume that equal KK number violation is possible in all internal directions, as in Eq. (25). In Appendix D, we show that for $n=2$ the number of decay channels for a KK mode of mass $m$ is approximately

$$\left(\frac{m_{h}}{m_{\rm KK}}\right)^{3/2}\delta n^{5/2}~.$$ (27)

Consequently, the total decay rate to KK gravitons is

$$\Gamma^{\rm tot}(m)\simeq\beta^{2}\frac{m_{h}^{4}\delta n^{3}}{M_{\rm Pl}^{2}m_{\rm KK}}~,$$ (28)

which corresponds to a typical lifetime of

$$\frac{1}{\Gamma^{\rm tot}}\simeq 0.1\,{\rm yr}\left(\frac{10}{\beta}\right)^{2}\left(\frac{10}{\delta n}\right)^{3}\left(\frac{10\,{\rm MeV}}{m_{h}}\right)^{4}\left(\frac{m_{\rm KK}}{{\rm eV}}\right)~.$$ (29)

To calculate precise limits and predictions, one must account for the distribution of KK graviton energies and masses, as well as time-dilation, but for our present purposes a rougher analysis is sufficient. For $n\geq 2$ most of the KK gravitons produced have masses not too much smaller than their energies, because the number of KK modes increases with mass, so we neglect the effect of time-dilation. We also approximate the mass distribution of KK gravitons from a supernova as a delta-function peaked at $m_{h}(t)$. Moreover, we assume that each KK graviton decays into two products of equal mass, each of mass $m_{h}/2$ (as shown in Appendix D, this is a reasonable approximation for the majority of decays). The typical graviton mass as a function of time then satisfies

$$\frac{dm_{h}(t)}{dt}\simeq-\Gamma^{\rm tot}(m_{h})\frac{m_{h}(t)}{2}~.$$ (30)

Hence, given an initial population of KK gravitons of mass $m_{0}$, we have

$$m_{h}(t)\sim m_{0}\left(1+\Gamma^{\rm tot}(m_{0})t/2\right)^{-1/4}~,$$ (31)

and the corresponding comoving number density of KK gravitons is $n(t)\simeq n_{0}m_{0}/m_{h}(t)$, where $n_{0}$ is the initial number density.

Suppose that an observational signal arises from energy injection into the SM at a time $t_{\rm obs}$ after a supernova. Then the corresponding rate of energy transfer to the SM, $\dot{E}_{\rm SM}$, is

$$\dot{E}_{\rm SM}(t_{\rm obs}) \simeq n(t_{\rm obs})m_{h}(t_{\rm obs})\frac{\lambda^{2}m_{h}(t_{\rm obs})^{3}}{M_{\rm Pl}^{2}}~.$$ (32)

Relative to the rate in the absence of KK number violation, $\dot{E}_{0}$, we have

$$\frac{\dot{E}_{\rm SM}(t_{\rm obs})}{\dot{E}_{0}} \simeq\left(1+\Gamma^{\rm tot}(m_{0})t_{\rm obs}/2\right)^{-3/4}~.$$ (33)

The typical timescale relevant to the heating of old neutron stars is roughly $10^{5}$ years, similar to that for direct searches for the decay products of gravitationally bound KK gravitons. The corresponding timescale for all past supernovae is even longer. Since these observations provide limits on $m_{\rm KK}$ that are not dramatically stronger than those from SN 1987A cooling even in the case of decays only to the SM, from Eqs. (29) and (33) it follows that for $\delta n$ of at least a few these constraints are much weaker than cooling bounds.

An interesting exception is the signal from Cas A, which has the shortest relevant timescale $t_{\rm obs}\simeq 300\,{\rm yr}$. Although previous analyses give constraints only about as strong as those from SN 1987A cooling, this could be improved by future observations. From Eqs. (29) and (33), for $\delta n\simeq 10$ decays to other KK gravitons render this signal irrelevant. However, if $\delta n$ is only slightly larger than $1$ it could still be important.

Additionally, previous limits considered only KK gravitons produced by bremsstrahlung. However, the pion production process is likely to be particularly important for these searches, since (as shown in Figure 5) it yields a substantially harder energy spectrum than bremsstrahlung. Consequently, the typical KK graviton mass is larger by roughly a factor of three, leading to decay rates (both intra-tower and to SM particles) that are more than an order of magnitude faster. The resulting SM decay products will also be more energetic, and therefore easier to detect above the background. Together, these effects could significantly enhance the sensitivity of future searches.

Meanwhile, for $n>2$ the number of available decay channels for a given KK graviton mass is even larger, and the signals from decays to the SM are correspondingly weaker than those from cooling for $\delta n$ non-zero.

We also note that an alternative realisation of KK number violation has been proposed and studied in detail in Ref. [Mohapatra:2003ah], in which the KK modes of a bulk scalar acquire expectation values. This leads to a different pattern of allowed decays than we have assumed. It might be interesting to explore the implications of this setup for the dark dimension scenario.

7 Discussion

In this paper we have revisited astrophysical cooling constraints on KK gravitons, motivated by recent interest in the dark dimension scenario. After analysing a wide range of astrophysical systems and considering two novel production mechanisms, we conclude that the bound from SN 1987A remains the most stringent. In this environment, even with conservative assumptions, the rate of production of KK gravitons from pions exceeds that from bremsstrahlung by a factor of a few, provided the soft approximation is accurate for the latter. As can be seen from Table 5, uncertainty in the protoneutron star profile of SN 1987A is a critical factor. Because of several partial cancellations, resonant production is relatively unimportant in supernovae, and it leads to limits from red giants that are a few orders of magnitude weaker than those from SN 1987A. Cold neutron stars also lead to constraints that are somewhat weaker than SN 1987A for $2$ and $3$ extra dimensions.

There are several ways that the robustness of the SN 1987A limit could be improved, although these are likely to be technically involved. The calculation of the bremsstrahlung rate using the soft approximation could be refined by analysing the impact of multi-body interactions or evaluating corrections that are formally higher order in the expansion but might be numerically significant. The rate could also be computed by evaluating the matrix elements directly, modifying the one-pion-exchange approximation with corrections or more sophisticated scattering potentials, building on work done in the context of axions, e.g., in Ref. [Carenza:2019pxu]. It would be especially useful to compute corrections to the leading-order chiral perturbation theory prediction for the pion process. On the simulation side, given the importance of pions for both KK graviton and axion limits, the inclusion of these in simulations would be valuable. Ultimately, full three-dimensional simulations including new sources of energy loss would be ideal, as in e.g. [Betranhandy:2022bvr]. However, significant progress is required even in the purely SM case, and the relevant BSM microphysics must be reliably determined before such extensions are worthwhile. The production of KK gravitons might also lead to observable effects in the evolution of supernova, beyond conventional cooling [Caputo:2025aac].

We have also shown that if decays within the KK graviton tower proceed at the rate typically assumed in the dark dimension scenario, searches for SM decay products lead to limits that are weaker than those from cooling bounds (although they remain significant if KK number violation is highly suppressed). Among the various proposed signals, the most promising are searches for SM particles from around very young neutron stars, where the population of heavy KK gravitons has had the least time to deplete. This is encouraging, since future observations of Cas A might provide sensitivity to previously unconstrained theories (the remnant of SN 1987A could also be observed, although it is further from Earth). As discussed, the pion production process could be especially important in this context, and a detailed study would be interesting.

For one extra dimension, cooling constraints from old neutron stars are comparable to those from SN 1987A, but both are substantially weaker than laboratory bounds. In this case, the best prospects for stronger future limits likely come from improved experiments searching for deviations from Newtonian gravity (on the other hand, one can hope for a discovery in specific non-minimal models that yield stronger signals, e.g. if there are neutrinos [Anchordoqui:2024xvl, Elacmaz:2025ihm, Eller:2025lsh] or axions [Mario] propagating in the 5 dimensional bulk). It is interesting to note that improved short-range gravity experiments could strengthen the lower bound on $m_{\rm KK}$ by about an order of magnitude, which would probe part of the dark dimension parameter space (see, e.g., [Manley:2024svv]). Meanwhile, for two extra dimensions, a possible future $100\,{\rm TeV}$ particle collider might be sensitive to theories currently allowed by SN 1987A constraints.

A notable omission from our present work is constraints from white dwarfs. Such systems lead to relatively strong bounds on axions, both from the luminosity function [Raffelt:1985nj, Blinnikov:1994eoa, MillerBertolami:2014rka, Isern_2018] and from the observed periods of pulsating white dwarfs [Corsico:2012sh, Corsico:2016okh, Battich:2016htm, Isern:2010wz]. Limits on KK gravitons have been obtained from pulsating white dwarfs in the past [Biesiada:2001iy], and it would be interesting to revisit these, although the previous limits were weaker than those from the cooling of SN 1987A. The resonant production process that we have considered might play a role, especially at high luminosities where neutrino cooling proceeds mainly due to plasmon decay. A reliable limit from the luminosity function would require implementing KK graviton emission in a white dwarf evolution model to capture the effects of feedback [MillerBertolami:2014rka], as well as careful analysis of the microphysics (e.g. ionic effects on bremsstrahlung emission [Nakagawa:1987pga]). Given that observational surveys are providing new data on the luminosity function [2015yCat..74464078K], this would be worthwhile. Pulsating white dwarfs are especially promising for KK gravitons since they can be relatively hot, and again new data is expected. It has also recently been shown that magnetic fields, such as those around magnetic white dwarfs, can affect resonant conversion of axions [Caputo:2020quz, Caputo:2021kcv, Brahma:2025vdr], and it might be interesting to investigate the consequences of this for KK gravitons.

Finally, although not the focus of our present work, several aspects of the cosmology of KK gravitons merit a brief mention. If KK gravitons produced from the SM thermal bath in the early Universe are effectively stable, then for $n=2$, a low reheating temperature of order $5\,{\rm MeV}$ is required to avoid overclosing the Universe, close to the minimum that is consistent with Big Bang Nucleosynthesis. If there are decays within the KK graviton tower with small $\delta n$, the decay products are approximately non-relativistic, and this bound on $T_{\rm RH}$ is essentially unchanged. For larger $\delta n$, the decay products are too warm to constitute all of dark matter. However, if $\delta n$ is sufficiently large, the decays might be fast enough that the KK graviton population remains relativistic after production, causing its energy density to redshift like radiation. It would be interesting to explore whether such scenarios can lead to viable cosmological histories with the KK gravitons as a warm, subdominant, component of dark matter. This population may, for example, contribute to $\Delta N_{\rm eff}$ at a level detectable with the forthcoming CMB-S4 experiment. It would also be interesting to study the signals of decays of dark matter KK gravitons to the SM in the $n=2$ case, analogous to the $n=1$ analysis of Ref. [Law-Smith:2023czn].

Appendix A Photon–KK graviton mixing

The photon–KK graviton mixing receives contributions from each species of charged fermion in the plasma $\Pi^{\mu,\alpha\beta}_{\gamma h}(Q)=\sum_{f}\Pi^{\mu,\alpha\beta}_{\gamma h,f}(Q)$. Using the Feynman rules in Ref. [Han:1998sg], and performing a calculation in the framework of real-time thermal field theory as described in Ref. [Hardy:2024gwy], we find that the mixing due to a fermion of charge $q_{f}$ with energy-momentum relation $E_{p}^{2}=p^{2}+m_{f}^{2}$ is given by

$$\Pi^{\mu,\alpha\beta}_{\gamma h,f}(Q)=q_{f}e\kappa\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{2E_{p}}\left(f_{e}(E_{p})-f_{\bar{e}}(E_{p})\right)\left(A^{\mu\alpha\beta}+\frac{B^{\mu\alpha\beta}+C^{\mu\alpha\beta}}{4(P\cdot Q)^{2}-Q^{4}}\right),$$ (34)

where

$$A^{\mu\alpha\beta}=\frac{2Q^{\mu}}{m_{A}^{2}\xi}(P^{\alpha}Q^{\beta}+P^{\beta}Q^{\alpha}-g^{\alpha\beta}(P\cdot Q))-(P^{\alpha}g^{\mu\beta}+P^{\beta}g^{\mu\alpha}),$$ (35)

$$B^{\mu\alpha\beta}=2 P^{\mu}\left[g^{\alpha\beta}\left(4(P\cdot Q)^{2}-Q^{4}\right)-4(P\cdot Q)\left(P^{\alpha}Q^{\beta}+P^{\beta}Q^{\alpha}\right)+Q^{2}Q^{\alpha}Q^{\beta}\right]$$ (36)

$$+ Q^{\mu}\left[Q^{2}\left(P^{\alpha}Q^{\beta}+P^{\beta}Q^{\alpha}\right)-8(P\cdot Q)P^{\alpha}P^{\beta}\right]$$ (37)

$$- Q^{2}(P\cdot Q)\left(Q^{\alpha}g^{\mu\beta}+Q^{\beta}g^{\mu\alpha}\right),$$ (38)

and

$$C^{\mu\alpha\beta}=4(P\cdot Q)^{2}\left(P^{\alpha}g^{\mu\beta}+P^{\beta}g^{\mu\alpha}\right)+8Q^{2}P^{\mu}P^{\alpha}P^{\beta}.$$ (39)

Here we have introduced a photon mass $m_{A}$ which is necessary to regulate the $t$ channel diagram, and $\xi$ is the usual gauge fixing parameter appearing in the photon’s propagator. The photon’s Lorentz index is $\mu$, and the KK graviton’s Lorentz indices are $\alpha\beta$. The $A^{\mu\alpha\beta}$ term comes from the $t$ channel and four-point contact diagrams shown in Table 1, while the $B^{\mu\alpha\beta}$ and $C^{\mu\alpha\beta}$ terms both come from the s-channel and u-channel diagrams.

Eq. (34) will be contracted with physical photon and KK graviton polarization vectors, which have the properties

$$\epsilon^{\gamma}_{\mu}Q^{\mu}=0,\qquad\epsilon^{h}_{\alpha\beta}Q^{\alpha}=\epsilon^{h}_{\alpha\beta}Q^{\beta}=0,\qquad\epsilon^{h}_{\alpha\beta}g^{\alpha\beta}=0.$$ (40)

Therefore it is easy to see that the $B^{\mu\alpha\beta}$ term will not contribute to the mixing. Furthermore, the $A^{\mu\alpha\beta}$ term reduces to $-(P^{\alpha}g^{\mu\beta}+P^{\beta}g^{\mu\alpha})$ so $\xi$ will not enter in any physical quantities, and the remaining terms are proportional to

$$q_{f}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{P^{\alpha}}{2E_{p}}\left(f_{e}(E_{p})-f_{\bar{e}}(E_{p})\right).$$ (41)

By rotational symmetry, only the $\alpha=0$ component can survive, however $P^{0}=E_{p}$ which cancels the term in the denominator and so the integral is simply proportional to $q_{f}(n_{f}-n_{\bar{f}})$. When we sum over the different fermion species, this gives the total charge density in the plasma, which must be zero as we have assumed that the plasma is in equilibrium. This means that $A^{\mu\alpha\beta}$ also does not contribute. Only the s-channel and u-channel contribute to the mixing, and only the $C^{\mu\alpha\beta}$ term remains. We can see that this is due to the presence of a tadpole in the thermal loop diagrams corresponding to the other channels.

Before we move on to perform the integration over momenta, we note some features of $\Pi^{\mu,\alpha\beta}_{\gamma h}(Q)$. Firstly we have the following Ward identities, also identified in Ref. [Gill:2023kyz], corresponding to the electromagnetic U(1) symmetry and gravitational diffeomorphisms

$$Q_{\mu}\epsilon^{h}_{\alpha\beta}\Pi^{\mu,\alpha\beta}_{\gamma h}(Q)=0$$ (42)

$$\epsilon^{\gamma}_{\mu}Q_{\alpha}Q_{\beta}\Pi^{\mu,\alpha\beta}_{\gamma h}(Q)=0.$$ (43)

Naively, contracting every index with the external momentum appears to give a non-zero value; however, once again after summing over all the fermion species in the plasma, charge neutrality enforces the identity

$$Q_{\mu}Q_{\alpha}Q_{\beta}\Pi^{\mu,\alpha\beta}_{\gamma h}(Q)=0.$$ (44)

These identities ensure that contributions at order $1/m_{h}$ vanish. While we expect this to hold to higher orders in perturbation theory, we have only checked the leading order case. It can also be checked that mixing only occurs between modes with the same helicity. For example the longitudinal polarization state of the photon does not mix with the helicity-1 states of the graviton, $\epsilon_{\mu}^{0}\epsilon_{\alpha\beta}^{\pm}\Pi^{\mu,\alpha\beta}_{\gamma h}=0$. To evaluate the integral, we assume $4(P\cdot Q)^{2}-Q^{4}\approx 4(P\cdot Q)^{2}$, originally argued for in Ref. [Braaten:1993jw]. This means that the denominator cancels in the first term of $C^{\mu\alpha\beta}$, leaving something of the same form as Eq. (41), which vanishes. The only part which contributes to the final mixing is therefore

$$\Pi^{\mu,\alpha\beta}_{\gamma h,f}(Q)_{\text{reduced}}=2Q^{2}q_{f}e\kappa\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{2E_{p}}\left(f_{e}(E_{p})-f_{\bar{e}}(E_{p})\right)\frac{P^{\mu}P^{\alpha}P^{\beta}}{(P\cdot Q)^{2}}.$$ (45)

After contracting with appropriate polarization vectors, the angular integration can be performed explicitly, yielding expressions of the form

$$\Pi^{L,T}_{\gamma h,f}(Q)=q_{f}e\kappa\int_{0}^{\infty}p^{2}\left(f_{e}(E_{p})-f_{\bar{e}}(E_{p})\right)H^{L,T}(E_{p},p,\omega,k)dp,$$ (46)

where $Q=(\omega,\bm{k})$ and $L,T$ denotes the mixing for the longitudinal and transverse modes respectively. We are free to shift our definition of $H^{L,T}$ by any constant $H^{L,T}\rightarrow H^{L,T}+f(\omega,k)$ since this will change the mixing by an amount proportional to the charge density which cancels out upon summing over species. We find that a convenient choice is to subtract off the zero-momentum contribution and define

$$\bar{H}^{L,T}(E_{p},p,\omega,k)=H^{L,T}(E_{p},p,\omega,k)-H^{L,T}(m_{f},0,\omega,k).$$ (47)

This both removes the unphysical $1/m_{h}$ divergence in the contribution from each individual species, and makes it easier to keep track of which species are important, since now in the non-relativistic limit $\bar{H}^{L,T}\propto p^{2}/m_{f}^{2}+\mathcal{O}(p^{4}/m_{f}^{4})$, meaning that the non-relativistic contribution is proportional to the pressure, which is suppressed by $T/m_{f}$. Since the $1/m_{h}$ terms only cancel out when overall charge neutrality is assumed, it is interesting to ask if such terms can exist out of equilibrium. These would lead to a large enhancement of the KK graviton production rate, and therefore to interesting observational signatures. We expect that the answer is no, and that the cancellation is enforced by Ward identities which rely only on the gauge symmetries of the theory. For a detailed discussion of the Ward identities and associated power counting see Ref. [Chivukula:2023qrt].

A.1 Summary of mixing calculation

We now present the final results of the calculation of the photon–KK graviton mixing $\Pi_{\gamma h,f}^{L,T}$ in various regimes. If the fermion, $f$, is non-degenerate and non-relativistic so that $p/m_{f}\ll 1$ for typical momenta, we find

$$\Pi_{\gamma h,f}^{L}(\omega)=-\frac{1}{\sqrt{6}}q_{f}e\kappa\frac{n_{f}T}{m_{f}}\frac{m_{h}}{\omega}v(3v^{2}-2),$$ (48)

$$\Pi_{\gamma h,f}^{T}=-\frac{1}{\sqrt{2}}q_{f}e\kappa\frac{n_{f}T}{m_{f}}\frac{m_{h}}{\omega}v,$$ (49)

where we have simplified the result using $\omega^{2}-k^{2}=m_{h}^{2}$, and introduced the graviton velocity $k=v\omega$. For the case where the fermions are ultra-degenerate, so that the distribution function can be approximated as a step function at the Fermi momentum $p_{F}$, the mixing is given by

$$\Pi_{\gamma h,f}^{L}=\frac{q_{f}e\kappa}{6\sqrt{6}\pi^{2}}p_{F}^{3}\frac{m_{h}}{v\omega}\left((3-2v^{2})H(v^{2}v_{F}^{2})-v^{2}v_{F}^{2}(1+H(v^{2}v_{F}^{2})\right),$$ (50)

$$\Pi_{\gamma h,f}^{T}=-\frac{q_{f}e\kappa}{3\sqrt{2}\pi^{2}}p_{F}^{3}\frac{m_{h}}{v\omega}\frac{G(v^{2}v_{F}^{2})}{2},$$ (51)

where $H(x)$ and $G(x)$ are functions which commonly appear in photon self-energy calculations [Raffelt:1996wa, Braaten:1993jw] and are defined as

$$G(x)=\frac{3}{x}\left(1-\frac{2x}{3}-\frac{1-x}{2\sqrt{x}}\log\left({\frac{1+\sqrt{x}}{1-\sqrt{x}}}\right)\right),\qquad H(x)=\frac{G(x)-x}{x-1}.$$ (52)

The Fermi velocity $v_{F}$ is defined as $p_{F}=v_{F}E_{F}$. Provided that the species in question also gives the dominant contribution to the photon’s self energy, the logarithms appearing in $H(x)$ and $G(x)$ can be eliminated by using the resonance condition $\Pi^{L/T}_{\gamma\gamma}(\omega_{*},k_{*})=\omega_{*}^{2}-k_{*}^{2}=m_{h}^{2}$ to give

$$\Pi_{\gamma h,f}^{L}=\frac{q_{f}e\kappa}{6\sqrt{6}\pi^{2}}p_{F}^{3}\frac{\omega_{*}m_{h}}{v_{*}\omega_{p}^{2}}\left((3-2v_{*}^{2})\left(1-\frac{\omega_{p}^{2}}{\omega_{*}^{2}}\right)-v_{*}^{2}v_{F}^{2}\right),$$ (53)

$$\Pi_{\gamma h,f}^{T}=\frac{q_{f}e\kappa}{3\sqrt{2}\pi^{2}}p_{F}^{3}\frac{m_{h}}{v_{*}\omega_{*}}\left(1-\frac{m_{h}^{2}}{\omega_{p}^{2}}\right),$$ (54)

where $v_{*}=k_{*}/\omega_{*}$. In all cases studied here, the electrons dominate both the mixing and the photon self-energy. However, one may imagine a situation in which one species is dominant in the photon self-energy, and a different species is dominant in the mixing. In this case Eqs. (53) and (54) would not apply. Notice that if the dominant species in the plasma is non-relativistic, such that $v_{F}\approx 0$, the longitudinal resonance satisfies $\omega_{*}\approx\omega_{p}$ and the transverse resonance satisfies $m_{h}\approx\omega_{p}$. In this limit both expressions vanish, as expected to match the non-relativistic case at lowest order in $T/m_{f}$. Making no approximations about degeneracy, but taking the fermions to be ultra-relativistic, we find after eliminating the logarithms with the resonance condition (now taking $v_{0}=1$)

$$\Pi_{\gamma h,f}^{L}=\frac{q_{f}e\kappa n_{f}}{2\sqrt{6}}\frac{\omega_{*}m_{h}}{v_{*}\omega_{p}^{2}}\left((3-2v_{*}^{2})\left(1-\frac{\omega_{p}^{2}}{\omega_{*}^{2}}\right)-v_{*}^{2}\right),$$ (55)

$$\Pi_{\gamma h,f}^{T}=\frac{q_{f}e\kappa n_{f}}{\sqrt{2}}\frac{m_{h}}{v_{*}\omega_{*}}\left(1-\frac{m_{h}^{2}}{\omega_{p}^{2}}\right),$$ (56)

which agrees with taking the relativistic limit of the ultra-degenerate expression. We again assume that the fermion, $f$, is the dominant contributor to the photon’s self energy.

If all particles are ultra-relativistic, then all contributions to the mixing are proportional to $q_{f}n_{f}$ so the leading order result shown here once again cancels out due to charge neutrality. Higher order terms are then suppressed by a factor of $m_{f}^{2}/T^{2}$, which is small in the ultra-relativistic limit.

Appendix B Red giant profiles

The profiles for the cores of stars at the tip of the red giant branch (with typical radius $10^{-2}R_{\odot}$) that we make use of when calculating constraints in Section 3 are shown in Fig. 6.

(a)

(b)

(c)

Figure 6: The density $\rho$, temperature $T$, and plasma frequency $\omega_{p}$ as a function of radial coordinate $r$ in the red giant cores for the 3 different metallicities, $Z_{\rm ini}$, that we use.

Appendix C Matrix element for the pion process

Here we give expressions for the matrix elements for KK graviton production from pions, corresponding to the Feynman diagrams in Figure 4. With the conventions for the momenta in that figure, we have

$$M^{s}_{\alpha\beta} =\frac{i\kappa f_{pp}}{4\sqrt{2}m_{\pi}}\frac{1}{(p_{1}+p_{2})^{2}-m_{N}^{2}}\left[(p_{1}+p_{2}+p_{4})^{\sigma}(\gamma_{\alpha}\eta_{\beta\sigma}+\gamma_{\beta}\eta_{\alpha\sigma}-2\gamma_{\sigma}\eta_{\alpha\beta})+4m_{N}\eta_{\alpha\beta}\right]$$ (57)

$$\hskip 140.0pt\times(\not{p_{1}}+\not{p_{2}}+m_{N})\not{p_{1}}\gamma^{5}$$

$$M^{u}_{\alpha\beta} =\frac{i\kappa f_{pp}}{4\sqrt{2}m_{\pi}}\frac{1}{(p_{2}-p_{3})^{2}-m_{N}^{2}}\not{p_{1}}\gamma^{5}(\not{p_{2}}-\not{p_{3}}+m_{N})$$ (58)

$$\hskip 140.0pt\times\left[(2p_{2}-p_{3})^{\sigma}(\gamma_{\alpha}\eta_{\beta\sigma}+\gamma_{\beta}\eta_{\alpha\sigma}-2\gamma_{\sigma}\eta_{\alpha\beta})+4m_{N}\eta_{\alpha\beta}\right]$$

$$M^{t}_{\alpha\beta} =\frac{i\kappa f_{pp}}{\sqrt{2}m_{\pi}}\frac{1}{(p_{1}-p_{3})^{2}-m_{\pi}^{2}}(\not{p_{1}}-\not{p_{3}})\gamma^{5}\left[p_{1}^{\rho}(p_{1}-p_{3})^{\sigma}C_{\alpha\beta\rho\sigma}+m_{\pi}^{2}\eta_{\alpha\beta}\right]$$ (59)

$$M^{c}_{\alpha\beta} =\frac{-i\kappa f_{pp}}{2\sqrt{2}m_{\pi}}p_{1}^{\sigma}(\gamma_{\alpha}\eta_{\beta\sigma}+\gamma_{\beta}\eta_{\alpha\sigma}-2\gamma_{\sigma}\eta_{\alpha\beta})\gamma^{5},$$ (60)

where $C_{\mu\nu\rho\sigma}$ is defined in Ref. [Han:1998sg]. Here $f_{pp}=g_{A}m_{\pi}/(2F_{\pi})\approx 1$ with $g_{A}\approx 1.3$ the axial-vector coupling constant of pions and nucleons and $F_{\pi}\approx 93{\,{\rm MeV}}$ is the pion decay constant.

Appendix D Decay channels

In this Appendix we count the number of channels by which a KK graviton of mass $m_{h}\gg m_{\rm KK}$ can decay to two lighter KK gravitons, assuming that the total KK number can change by up to $\delta n$ during the decay, as described in Section 6.

Consider first the case of one extra dimension. We label the KK number of an initial state $N\equiv m_{h}/m_{\rm KK}$ (we take this to be positive without loss of generality), and those of the final state, $a$ and $b$. The number of decay channels is given by the number of integers $a$, $b$ such that

$$N-\delta n<a+b<N~,$$ (61)

where the lower bound is from the maximum violation of KK number, and the upper bound is from conservation of energy. Hence, the number of channels is of order

$$N\delta n=\frac{m_{h}}{m_{\rm KK}}\delta n~,$$ (62)

(we are not concerned with order-one numerical factors).

Next we consider the case of two equal radius extra dimensions, assuming violation of up to $\delta n$ is allowed for the KK numbers associated with both compactified directions. Without loss of generality, we take the KK numbers of the initial state to be $(N,0)$, where $N=m_{h}/m_{\rm KK}$ as before. Call the KK numbers of the final states ${\bf a}=(a_{1},a_{2})$ and ${\bf b}=(b_{1},b_{2})$. We must count the number of pairs of vectors over the integers $({\bf a},{\bf b})$ that satisfy

$$N-\delta n<a_{1}+b_{1}<N+\delta n\,,\quad{\rm and}\quad-\delta n<a_{2}+b_{2}<\delta n~,$$ (63)

(from KK number) and

$$|{\bf a}|+|{\bf b}|<N~,$$ (64)

from energy conservation. These conditions are equivalent to

$$|{\bf a}|+|{\bf f}-{\bf a}|<N~,$$ (65)

where ${\bf f}\equiv(f_{1},f_{2})={\bf a}+{\bf b}$ is such that $f_{1}\in(N-\delta n,N+\delta n)$ and $f_{2}\in(-\delta n,\delta n)$.

Consider a particular fixed ${\bf f}=(N-\delta n_{1},\delta n_{2})$. Eq. (65) corresponds to points inside an ellipse with foci at ${\bf 0}$ and ${\bf f}$ and major axis of length $N$. The number of such points is the area of the ellipse, which (dropping numerical factors) is

$$N\sqrt{N^{2}-|{\bf f}|^{2}}\simeq N\sqrt{N\delta n_{1}}~,$$ (66)

where we used $\delta n_{1}\leq\delta n\ll N$.

Since the number of different ${\bf f}$ is of order $\delta n^{2}$, the total number of open decay channels is

$$\delta n^{2}N\sqrt{N\delta n}=\left(\frac{m_{h}}{m_{\rm KK}}\right)^{3/2}\delta n^{5/2}~,$$ (67)

where we neglect an overall numerical factor that is about $0.3$. Eq. (67) is the expression used in the main text.