# Dark matter spikes in the vicinity of Kerr black holes

Francesc Ferrer, Augusto Medeiros da Rosa, Clifford M. Will

► To cite this version:

Francesc Ferrer, Augusto Medeiros da Rosa, Clifford M. Will. Dark matter spikes in the vicinity of Kerr black holes. Physical Review D, 2017, 96 (8), pp.083014. 10.1103/PhysRevD.96.083014. hal-01903337

HAL Id: hal-01903337

https://hal.science/hal-01903337v1

Submitted on 4 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dark matter spikes in the vicinity of Kerr black holes

Francesc Ferrer,¹ Augusto Medeiros da Rosa,¹ and Clifford M. Will²,³

¹Department of Physics, McDonnell Center for the Space Sciences, Washington University,
St. Louis, Missouri 63130, USA

²Department of Physics, University of Florida, Gainesville, Florida 32611, USA

³GReCO, Institut d'Astrophysique de Paris, CNRS, Université Pierre et Marie Curie,
98 bis Boulevard Arago, 75014 Paris, France

(Received 29 July 2017; published 20 October 2017)

The growth of a massive black hole will steepen the cold dark matter density at the center of a galaxy into a dense spike, enhancing the prospects for indirect detection. We study the impact of black hole spin on the density profile using the exact Kerr geometry of the black hole in a fully relativistic adiabatic growth framework. We find that, despite the transfer of angular momentum from the hole to the halo, rotation increases significantly the dark matter density close to the black hole. The gravitational effects are still dominated by the black hole within its influence radius, but the larger dark matter annihilation fluxes might be relevant for indirect detection estimates.

DOI: 10.1103/PhysRevD.96.083014

## I. INTRODUCTION

The central regions of galaxies and clusters of galaxies are a prime target for searches for indirect signals of annihilations or decays of dark matter (DM) particles. Arguably, the brightest DM source in the sky is the center of our Galaxy because of its close proximity and potentially large concentration of DM. Although the precise DM distribution is not well constrained by observations, measurements of the stellar rotation curve provide robust evidence for its presence in the innermost sector of our Galaxy [1]. Dissipationless DM-only simulations of halos with masses ranging from dwarf galaxies to rich clusters suggest that its density follows a near-universal cusped profile [2], which would enhance fluxes of high-energy radiation originated by DM reactions. Baryonic effects could modify the inner shape of the DM density profile, either steepening the profile via adiabatic contraction [3,4] or softening the cusp into a core through repeated and violent oscillations in the central potential due to energy injection from active galactic nuclei or supernovæ [5–8]. Although neither hydrodynamic simulations (see e.g. [9] for a review) nor dynamical constraints [10,11] can determine the exact DM brightness of the Galactic Center, several intriguing observations of possible signals have fuelled a sustained interest in understanding the DM distribution in the central regions of the Galaxy and of large scale structures in general [12–18].

In addition, there is strong evidence that the Galaxy harbors a massive black hole ($M_{BH} \gtrsim 4 \times 10^6 M_{\odot}$) at the center (see e.g. [19] for a review), which could lead to a significant increase of the DM density in its neighborhood. As the black hole grows and pulls in more dark matter, the density distribution becomes steeper. A Newtonian analysis with an *ad hoc* treatment of particle capture by the hole showed that a *spike* in the dark-matter density is created,

which causes a significant boost in the DM annihilation fluxes [20].

A fully relativistic calculation for the case of a spherical hole was completed in [21] (hereafter referred to as SFW), concluding that the Newtonian framework underestimates the dark matter density very close to the black hole: the spike reaches significantly higher densities, and it extends closer to the event horizon. Using the Schwarzschild geometry, SFW found a closed form for the boundary of the region in phase space containing bound orbits that do not cross the event horizon, and consistently took particle capture by the black hole into account.

On the other hand, a typical black hole is expected to be rotating rather fast. The spin of a super-massive hole depends on whether it gained most of its mass via mergers or accretion. BHs that grow mostly through disk accretion, adding material with constant angular momentum axis, end up spinning rapidly [22]. Although the outcome of an individual merger event depends on the initial spin alignment, for BHs that grow through repeated mergers we expect a distribution of spins that peaks at $\tilde{a} \sim 0.7$, where *a* is the Kerr parameter related to the angular momentum *J* by $a \equiv J/m$ [23–25] and $\tilde{a} \equiv a/Gm$ is the associated dimensionless quantity (we use units in which the speed of light $c = 1$). Millimeter VLBI observations of Sgr A* [26,27], and the analysis of Quasi-Periodic Oscillations of hot plasma spots in the surrounding orbiting material [28] suggest a value of $\tilde{a} \sim 0.65$ for the spin of the super-massive black hole in the Galactic Center, which obtains its angular momentum through accretion of tidally disrupted stars and gas clouds with randomly oriented angular momenta.

As pointed out in SFW, even if the initial DM distribution is spherically symmetric, the dragging of inertial frames induced by the rotation of the black hole could create a flux of DM in the azimuthal direction proportional
to the Kerr parameter $a$. Numerical investigations of dark matter geodesics around a preexisting hole have found interesting features in the annihilation spectrum [29]. We here study the modifications in the DM spike that ensue from the rotation of the black hole by extending the calculation in SFW to the Kerr geometry.

The main difficulty in performing this calculation for a spinning black hole is that, in general, no analytic expression is known for the critical Carter constant that separates orbits that plunge into the event horizon from those that do not. We circumvent this difficulty by using a brute-force determination of orbit “stability” from its turning points. We were able to find a closed form of the phase-space boundary for the subset of orbits that are contained in the equatorial plane, and we check our numerical calculations in this particular case.

The main effect we observe is a further enhancement of the density profile in its innermost region, at around 5 times the gravitational radius $Gm$ of the black hole. Our results are obtained from an initially symmetric dark matter distribution and can be understood as coming from the preferential binding of corotating orbits to the black hole. These particles “feel” a deeper potential well than the Schwarzschild case and are thus pulled closer to the hole, in what turns out to be the dominant effect over the preferential capture of counter-rotating particles by the black hole itself.

The increase obtained is more pronounced in the equatorial plane and is more relevant for high black hole spins. When $\tilde{a} = 0.8$, we obtain a peak density that is around 70% higher than the spherical case for DM that initially follows a Hernquist profile, as shown in Fig. 10. For a near-extreme black hole, the density just outside the ergosphere is more than an order of magnitude greater than the peak density obtained in the Schwarzschild case.

We choose the Hernquist profile as a proxy for a cuspy dark matter distribution as suggested by dark matter only N-body simulations, which do not include the baryonic component of the Universe. As suggested by hydrodynamical simulations [8,30–33], baryons could play an important role in determining the shape of the halo, possibly leading to the formation of a core. The results in Fig. 7 for the constant distribution provide a good description in this case, and show a 20% peak density increase in the equatorial plane for $\tilde{a} = 0.8$ compared to Schwarzschild.

Let us stress that we are assuming that the growth of the black hole is adiabatic, but several dynamical effects could affect our conclusions. For instance, if the black hole spirals from an initially off-center location [34], if significant merger events occur [35], or if gravitational scattering off stars heats the dark matter [36–38], the density inside the spike could be considerably lower. Recent observations for the case of Sgr A* show that the density of old stars is flat [39,40], or even decreasing towards the Galactic Center, which implies heating timescales well above 10 Gyr. However, it has been pointed out that this result only

holds for a small fraction of bright stars and the evidence for the existence of a central cusp in the vicinity of the Milky Way’s central black hole is at present inconclusive (see e.g. [41] for a recent review). In addition, as discussed below, the effects of dark matter annihilations could deplete and weaken the density profile [42,43]. Although these effects are important, our main purpose is to understand the general relativistic effects due to the rotation of the black hole, extending the nonrelativistic treatment in [20] and the relativistic static calculation in SFW.

On the other hand, it is important to note that observables such as fluxes depend on integrals of the density profile, and the region where the enhancement occurs has a very small volume. Thus, the impact of the enhancement on integrated effects will be small, but should still be taken into account in model building [14,44,45]. Moreover, a significant number of super-massive black holes in Active Galactic Nuclei are known to be rapidly spinning [46], which could potentially enhance their contribution to the isotropic gamma-ray background [47].

The rest of this paper provides the details supporting these conclusions. In Sec. II we describe the phase space available to bound orbits around a Kerr black hole. In Sec. III we obtain the DM density profile resulting from an initially constant phase-space distribution, and in Sec. IV we obtain the profile from an initially cuspy Hernquist distribution. Section V discusses the implications of these DM distributions for the gravitational environment around the black hole and for the fluxes of radiation from the DM spike. Concluding remarks are presented in Sec. VI.

## II. BLACK HOLE GROWTH IN A DARK MATTER HALO

We begin by reviewing the analysis of how the adiabatic growth of a black hole modifies the density and velocity dispersions within a preexisting dark matter halo. We follow the general relativistic approach in SFW, which extended the Newtonian treatment in [48–50] (see also [51]).

Our starting point is the relativistic phase space distribution $f^{(4)}(x, p)$ describing a system of dark matter particles of rest mass $\mu$ [52–54], normalized so that when integrated over phase space it gives the total mass of the halo. The mass current four-vector can be written in terms of the distribution function as

$$J^{\mu}(x) = \int f^{(4)}(x, p) u^{\mu} \sqrt{-g} d^4 p, \quad (1)$$

where $u^{\mu} = p^{\mu}/\mu$ is the four-velocity, and $g \equiv \det(g_{\mu\nu}(x))$ is the determinant of the metric.<fnref n="1" /> Knowledge of the mass

<footnote n="1">¹Note that we use contravariant components $dp^{\mu}$ to define the volume element in momentum space, but the same results can be obtained by taking covariant components $p_{\mu}$ as the argument of the distribution function.</footnote>
current four-vector will allow us to find the dark matter density: from the definition $J^\mu = \rho u^\mu$ and the fact that $u^\mu u_\mu = -1$, we obtain the mass density as measured in a local freely falling frame as $\rho = \sqrt{-J_\mu J^\mu}$.

To calculate $J^\mu$ we need to know the distribution function $f^{(4)}(x, p)$ and the boundary of integration over momentum space. Both tasks are greatly simplified if instead of $p^\mu$ we use invariant constants of the motion to write down the distribution function and the volume element. We review below this change of variables for the case of a Kerr black hole background following the discussion in SFW.

We use Boyer-Lindquist coordinates to express the Kerr line element for a hole of mass $m$ (with $c = 1$):

$$ds^2 = -\left(1 - \frac{2Gmr}{\Sigma^2}\right)dt^2 + \frac{\Sigma^2}{\Delta}dr^2 + \Sigma^2 d\theta^2 - \frac{4Gmar}{\Sigma^2}\sin^2\theta d\phi dt + \left(r^2 + a^2 + \frac{2Gmra^2\sin^2\theta}{\Sigma^2}\right)\sin^2\theta d\phi^2, \quad (2)$$

where $a$ is the Kerr parameter defined above, and we have introduced the functions $\Delta = r^2 + a^2 - 2Gmr$ and $\Sigma^2 = r^2 + a^2 \cos^2\theta$.

The orbits of dark matter particles of rest mass $\mu$ in this geometry admit four constants of the motion: the energy per unit mass, $\mathcal{E}$, and angular momentum per unit mass, $L_z$, that are conserved because the metric is stationary and axisymmetric; the mass-shell condition; and the so-called Carter constant per unit (mass)$^2$, $C$ [55]:

$$\begin{aligned} \mathcal{E} &\equiv -u_t = -g_{tt}u^t - g_{t\phi}u^\phi, \\ L_z &\equiv u_\phi = g_{\phi\phi}u^\phi + g_{t\phi}u^t, \\ \mu &\equiv \sqrt{-p_\mu p^\mu}, \\ C &\equiv \Sigma^4(u^\theta)^2 + \frac{L_z^2}{\sin^2\theta} + a^2\cos^2\theta(1 - \mathcal{E}^2). \end{aligned} \quad (3)$$

We use a definition of the Carter constant that, in the spherically symmetric Schwarzschild limit ($a \to 0$), reduces to the square of the conserved total angular momentum per unit mass, $L: C \to L^2 = u_\theta^2 + u_\phi^2/\sin^2\theta$.

Given these definitions and using $p^\mu = \mu u^\mu$, we calculate the (inverse) Jacobian

$$\mathcal{J} = \left| \frac{\partial(\mathcal{E}, C, L_z, \mu)}{\partial(p^t, p^r, p^\theta, p^\phi)} \right| = \frac{2\Sigma^4 \Delta |u_r||u^\theta|\sin^2\theta}{\mu^3}. \quad (4)$$

The last step is to write the necessary four-velocity components appearing in the Jacobian in terms of the constants of the motion. These are

$$\begin{aligned} u^\theta &= \pm \frac{1}{\Sigma^2} \sqrt{U(\theta)} \\ u_r &= \pm \frac{r^2}{\Delta} \sqrt{V(r)}, \end{aligned} \quad (5)$$

where

$$U(\theta) \equiv C - \frac{L_z^2}{\sin^2\theta} - a^2(1 - \mathcal{E}^2)\cos^2\theta, \quad (6)$$

and

$$V(r) \equiv \left(1 + \frac{a^2}{r^2} + \frac{2Gma^2}{r^3}\right) \mathcal{E}^2 - \frac{\Delta}{r^2} \left(1 + \frac{C}{r^2}\right) + \frac{a^2 L_z^2}{r^4} - \frac{4Gma\mathcal{E}L_z}{r^3}. \quad (7)$$

The presence of the $\pm$ signs in Eq. (5) implies that the $r$ and $\theta$ components of the current $J_\mu$ vanish, as they come with an absolute value in the Jacobian and we must integrate over both positive and negative values of $u_r$ and $u_\theta$ in Eq. (1), since they are equally likely to be positive or negative for a given set of values for $\mathcal{E}$, $C$ and $L_z$. For the other two components, we must put in an extra factor of 4 to allow for the integration over the positive and negative values of both four-velocity components. Note that, as pointed out in SFW, $J^\phi$ will not vanish even for spherically symmetric dark matter distributions, $f(\mathcal{E}, C)$, since the last term in Eq. (7) is linear in $L_z$. This effect, related to the dragging of inertial frames caused by the rotation of the black hole, will influence the dark matter spike around the black hole as we further elaborate below.

Concerning the distribution function, we are assuming that all of the dark matter particles have the same mass $\mu$. Then we are allowed to write $f^{(4)}(x, p) = \mu'^{-3}f(\mathcal{E}, C, L_z)\delta(\mu' - \mu)$, which will allow us to perform the integration over $d\mu'$.

Putting everything together, we can write the nonzero components of the mass current four-vector as

$$J_t(r, \theta) = \frac{-2}{r^2 \sin \theta} \int d\mathcal{E} dC dL_z \frac{\mathcal{E} f(\mathcal{E}, C, L_z)}{\sqrt{V(r)} \sqrt{U(\theta)}}, \quad (8a)$$

$$J_\phi(r, \theta) = \frac{2}{r^2 \sin \theta} \int d\mathcal{E} dC dL_z \frac{L_z f(\mathcal{E}, C, L_z)}{\sqrt{V(r)} \sqrt{U(\theta)}}, \quad (8b)$$

where we used $u_t = -\mathcal{E}$ and $u_\phi = L_z$ to show the covariant components.

Introducing $\Omega \equiv J_\phi/J_t$, the density in the rest frame of the distribution, $\rho = \sqrt{-g^{\mu\nu}J_\mu J_\nu}$, reads$^2$:

<footnote>$^2$We correct a typo in Eq. (3.11) in SFW.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. A plot of $V_{\text{eff}} \equiv -V(r)$ vs $r/Gm$ for a Kerr black hole with $\tilde{a} = 0.75$, $\mathcal{E} = 0.97$, $C = 12(Gm)^2$, and $L_z = 2Gm$. This is equivalent to the usual effective potential of classical mechanics. The two “bound orbits” are clearly seen ranging along the shaded regions. The bottom of the leftmost peak is not shown, as it is much lower than the normal bound orbit.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. A plot of $V_{\text{eff}}(r)$ vs $r/Gm$ with fixed $\mathcal{E}$, $L_z$. Varying the Carter constant separates bound and plunge orbits. The solid blue line depicts the same orbit shown in Fig. 1. Decreasing $C$ we reach the dashed magenta line showing the critical orbit with $C = C_{\text{crit}}$, which is equal to $C_{\text{crit}} = 10.3(Gm)^2$ in this case. For smaller $C$ the orbits plunge into the hole, as depicted by the dot-dashed red line with $C = 9(Gm)^2$.</figcaption></figure>

$$ \rho = |J_t| \sqrt{\frac{g_{\phi\phi} - 2g_{t\phi}\Omega + g_{tt}\Omega^2}{\Delta\sin^2\theta}}. \quad (9) $$

To actually evaluate Eq. (8), we still need to specify the shape of the integration region over $\mathcal{E}$, $C$ and $L_z$, as well as the distribution of the trajectories of dark matter particles $f(\mathcal{E}, C, L_z)$. The first task is quite involved and is one of the main results of this work. We give several general considerations here that are further developed in the following sections.

### A. Region of integration in $\mathcal{E}-C-L_z$ space

If we consider only the contribution of bound particles to the density, there are three constraints we must apply to our phase space: the energy $\mathcal{E}$ is bounded above by 1, and we must also have $V(r) \ge 0$ and $U(\theta) \ge 0$. Also, as pointed out in previous work [20,21], we must remove orbits that plunge into the black hole. For the case of a Schwarzschild black hole, this capture condition can be worked out analytically. The criterion obtained in SFW is that, given an energy, there is a critical value of the angular momentum $L_{\text{crit}}$ below which all orbits are captured.

For the Kerr metric, no such constraint has been derived. For given values of the conserved quantities $(\mathcal{E}, C, L_z)$, we now have in general two sets of “bound orbits”: one orbit with two turning points beyond the horizon, and one with a subhorizon turning point [56]. These can be seen extending through the shaded regions in Fig. 1, which shows the effective potential $V_{\text{eff}} \equiv -V(r)$. In between the two shaded regions corresponding to the two orbits, the effective potential has a maximum at $r_{\text{unst}}$.

We can thus exclude the plunge orbits on a case by case basis using the following criteria:

(i) The value of $r$ at which we are evaluating the current must be to the right of the unstable orbit $r_{\text{unst}}$.

(ii) The point $r_{\text{unst}}$ must also be in the forbidden region, $V_{\text{eff}}(r_{\text{unst}}) = -V(r_{\text{unst}}) > 0$, so that there is a potential barrier between the orbit of interest and the horizon. This will ensure that the orbit does not cross the horizon and become trapped.

In this way we exclude all plunge orbits and all unphysical orbits with a subhorizon turning point, as illustrated in Fig. 2. Unfortunately, it is not possible in general to find analytic expressions of the region in $(\mathcal{E}, C, L_z)$-space where these two criteria are met, so our code implements the capture condition numerically. Nevertheless, useful insight can be gained by focusing on a subset of the orbits. The case of nonrelativistic particles with $\mathcal{E} \approx 1$ was studied by Will [57], who found an approximate analytic expression for the critical value of $C$. We consider in Sec. III A all bound orbits that are contained in the equatorial plane. In this case $L_z = \pm\sqrt{C}$, which makes the problem tractable and also provides a useful check of the full numerical calculation.

### B. Adiabatic approximation

Now that the available phase space has been determined, at least implicitly, we turn to $f(\mathcal{E}, C, L_z)$. Our starting point is a known initial nonrelativistic dark matter distribution without a central black hole, $f'(E, L^2, L_z)$. For simplicity, we will assume our initial distribution to be spherically symmetric, generating a potential $\Phi(r)$.

As reviewed in appendix A, the adiabatic growth of the black hole preserves the form of the distribution function, $f(\mathcal{E}, C, L_z) = f'(E, L^2, L'_z)$. Here, $E, L^2$ and $L'_z$ are obtained from the Kerr constants of the motion by noting that each particle responds to the slow change in the gravitational potential by altering its energy $E$ and angular momentum $L$ and $L_z$, in such a way that the action variables $I_r, I_\theta$ and $I_\phi$ are kept fixed [48,58].
For an orbit in the initial nonrelativistic dark matter distribution, these are

$$ 
\begin{align} 
I'_r(E, L) &\equiv \oint dr \sqrt{2E - 2\Phi - L^2/r^2}, \\ 
I'_\theta(L, L'_z) &\equiv \oint d\theta \sqrt{L^2 - \frac{L'_z^2}{\sin^2\theta}} = 2\pi(L - |L'_z|), \\ 
I'_\phi(L'_z) &\equiv \oint d\phi L'_z = 2\pi L'_z. 
\end{align} 
\tag{10} $$

For a bound orbit in the Kerr geometry,

$$ 
\begin{align} 
I_r(\mathcal{E}, C, L_z) &\equiv \oint u_r dr = \oint dr \frac{\sqrt{V(r)}}{1 + \frac{a^2}{r^2} - \frac{2Gm}{r}}, \\ 
I_\theta(\mathcal{E}, C, L_z) &\equiv \oint u_\theta d\theta = \oint d\theta \sqrt{U(\theta)}, \\ 
I_\phi(L_z) &\equiv \oint u_\phi d\phi = 2\pi L_z. 
\end{align} 
\tag{11} $$

The equality $I_\phi = I'_\phi$ implies $L'_z = L_z$. The conservation of this component of the angular momentum should be expected as both the initial and the final states have axial symmetry.

Given a set $(\mathcal{E}, C, L_z)$ for a Kerr orbit, together with $L_z = L'_z$, we can determine $L$ via

$$ L = |L'_z| + \frac{I_\theta(\mathcal{E}, C, L_z)}{2\pi}. $$

Given this value of $L$, we can obtain $E$ by equating the radial actions. Since the radial integral in (10) cannot in general be solved analytically, we find the energy $E$ using the bisection method.

### III. DENSITY PROFILE FOR A CONSTANT DISTRIBUTION FUNCTION

The simplest possible example that can be considered is that of a constant distribution function. Although somewhat unrealistic, it can be seen as a toy model for describing the stars close to a black hole forming in the core of an approximately isothermal system. This case does not require any adiabatic matching, but is still useful for building intuition, as the current-density will be directly related to the total phase space volume.

#### A. Restriction to equatorial orbits

We first consider the subset of particles following planar orbits, which are only possible on the equatorial plane. In this case, the calculations can be carried out analytically to a large extent. Moreover, the result will be a lower bound on the final dark matter density on the plane, because the calculation omits nonplanar orbits that cross $\theta = \pi/2$.

We will begin by setting $f(\mathcal{E}, C, L_z) = f_{\text{eq}}\delta(u_\theta)$ and focusing on $\theta = \pi/2$. This allows us to perform the integrals in Eq. (8) over $C$ and $L_z$, obtaining

$$ 
\begin{align} 
J_t &= \frac{-8f_{\text{eq}}}{\sqrt{Gmr}} \int_{\mathcal{E}_{\text{min}}(r)}^{1} d\mathcal{E} \mathfrak{F}_t(\mathcal{E}), \\ 
J_\phi &= \frac{8f_{\text{eq}}\sqrt{Gm}}{\sqrt{r}} \int_{\mathcal{E}_{\text{min}}(r)}^{1} d\mathcal{E} \mathfrak{F}_\phi(\mathcal{E}). 
\end{align} 
\tag{12} $$

The explicit forms of the functions $\mathfrak{F}_{t,\phi}$ representing the $C-L_z$ integration are given in Eq. (B17). The $\mathcal{E}$ dependence in these functions is implicit in the quantities $L_z^\pm$, which are the two roots of the equation $V(r) = 0$ given in Eq. (B5), and $L_z^{\text{crit}}$, the critical angular momentum for capture by the black hole. For a corotating planar orbit, $L_z^+ \ge L_z^{\text{crit,+}} > 0$, while for a counter-rotating orbit $L_z^- \le L_z^{\text{crit,-}} < 0$. The currents are thus naturally separated in a counter-rotating and a corotating part corresponding to the integration domain $\Delta L_z$ in $L_z$:

$$ \Delta L_z = (L_z^-, L_z^{\text{crit,-}} \cup (L_z^{\text{crit,+}}, L_z^+)). $$

The minimum energy $\mathcal{E}_{\text{min}}(r)$ is found by setting $L_z^{\text{crit}\pm}(\mathcal{E}) = L_z^{\pm}(\mathcal{E})$.

As further discussed in appendix B, this is a useful simplification because the numerical difficulty of the calculation is transferred to the pair of functions $L_z^{\text{crit}\pm}(\mathcal{E})$. Results for the density and the co- and counter-rotating parts of $J_t$ are shown in Figs. 3 and 4. Figure 3 shows that increasing spin increases the density, and Fig. 4 is helpful in understanding the physical origin of this effect.

As Fig. 4 shows, the density enhancement is coming from the corotating orbits, which are more deeply bound to the black hole. Since this binding energy increases with the spin parameter, the density will also increase with black hole spin. This is illustrated in Fig. 5, which shows the

<figure><img src="image_4.png" /><figcaption>FIG. 3. Density profiles for a distribution of equatorial orbits. The density increases as we vary the Kerr parameter $\tilde{a} = 0$ (solid), 0.5 (dashed), 0.8 (dot-dashed) and 0.998 (dotted). Since the latter value is greater than $\tilde{a} > 2(\sqrt{2} - 1)$, the spike extends into the ergosphere (shaded region).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. The corotating (dashed) and counter-rotating (dotted) parts of $|J_t|$ for $\tilde{a} = 0.8$, compared with $|J_t|/2$ for the Schwarzschild case (solid).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 6. The blue region shows a phase space slice of fixed energy, $\mathcal{E} = 0.98$, at $r/Gm = 20$ on the equatorial plane for a Kerr hole with $\tilde{a} = 0.5$. The upper dashed line corresponds to $V(r) = 0$ for $r = 20Gm$, and the lower dashed line corresponds to the capture condition for $\mathcal{E} = 0.98$, both in the Schwarzschild case. The red-shaded lower-left region shows the counter-rotating orbits that are lost due to capture by the hole, which are compensated by the tightly bound corotating orbits in the dark blue-shaded region.</figcaption></figure>

minimum allowed energy for $r/Gm = 10$ as a function of spin parameter $\tilde{a}$.

As explained in SFW, the density is zero at the coordinate $r$ such that $\mathcal{E}_{\min}(r) = 1$. Orbits that go any closer to the black hole will have to be unbound in order to not be captured, and are therefore not included in this calculation. For the equatorial plane, this occurs at $r_{\text{mb}}^{\pm}/Gm = 2 \mp \tilde{a} + 2\sqrt{1 \mp \tilde{a}}$, which corresponds to the radius of the marginally bound circular orbit, where the upper (lower) sign is for corotating (counter-rotating) orbits. Note that the density vanishes at $r_{\text{mb}}^{\pm}$, which can be inside the ergosphere for $\tilde{a} \ge 2(\sqrt{2}-1) \approx 0.83$. Only corotating orbits contribute in the range $r_{\text{mb}}^{+} \le r \le r_{\text{mb}}^{-}$. In a small region around $r \gtrsim r_{\text{mb}}^{-}$, when counter-rotating orbits start contributing, the slope of the density profile becomes less steep. For Schwarzschild, $r_{\text{mb}}^{+} = r_{\text{mb}}^{-} = 4Gm$, which is where SFW found the end point of the spike in this case.

It is important to note that the boost obtained here is amplified by the fact that the corotating equatorial orbits are the most bound to the black hole. Nevertheless, when we integrate over the full phase space, the net effect is still the same: the enhanced binding of the corotating orbits is large enough to make up for the loss of counter-rotating orbits, which are preferentially captured by the black hole. This

last point is also illustrated in Fig. 6, which shows a fixed energy slice on the equatorial plane. The planar orbits fall on the left and right boundaries of the blue region, but the trade-off between loss of counter-rotating orbits that is more than compensated by the addition of corotating ones holds in general.

## B. Full phase space

Unlike the equivalent calculation for the Schwarzschild case done in SFW, there is no analytic expression for the integration volume in the Kerr geometry. None of the integrals can be performed analytically to simplify the current density. Therefore, we used a Monte Carlo (MC) method to overcome these difficulties. To maximize efficiency in the MC evaluation of the integrals, we look for the $(\mathcal{E}, C, L_z)$-cube that has the tightest fit to our phase space: we know that the maximum value of $\mathcal{E}$ is 1. We find an upper bound for $C$ by noting that, from the positivity of $U(\theta)$ one obtains $C \ge L_z^2$ for bound orbits.

Thus, substituting $L_z = -\sqrt{C}$ and $\mathcal{E} = 1$ in $V(r)$, one can obtain an upper bound on $V(r)$ and, consequently, on $C$, as long as $r/Gm > 2$. For $r/Gm < 2$ or, in general, for any $r$ within the ergoregion of the black hole, this upper bound on $C$ can be found by looking for a plunge orbit with $\mathcal{E} = 1$ and $V(r) = 0$.

Once we have the upper bound on $C$ and $\mathcal{E}$, we find the orbit with minimum energy and Carter constant $(\mathcal{E}_{\min}(r, \theta), C_{\text{crit}}(r, \theta))$ by requiring that all of the bounds in our phase space be satisfied simultaneously. This means $V(r) = U(\theta) = 0$, and this orbit must also be a plunge orbit, implying that the potential will have a double root at its unstable orbit $r_{\text{unst}}$. This gives us a system of four polynomial equations for $(\mathcal{E}_{\min}, C_{\text{crit}}, L_z^*, r_{\text{unst}})$, which we solve using a homotopy continuation method [59].

<figure><img src="image_7.png" /><figcaption>FIG. 5. Minimum allowed energy for $r/Gm = 10$ as a function of Kerr parameter $\tilde{a}$ in the equatorial plane.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 7. The dark matter density in the equatorial plane increases with the spin parameter, and the spike gets closer to the hole. The different lines show the Schwarzschild calculation from SFW (blue, dot-dashed), $\tilde{a} = 0.5$ (purple, dashed), and $\tilde{a} = 0.8$ (red, solid).</figcaption></figure>

We find it advantageous to implement the positivity of $U(\theta)$ explicitly through the change of variables:

$$ 
\begin{aligned}
\mathcal{E} &= x + (1-x)\mathcal{E}_{\min}, \\
C &= yC_{\max} + (1-y)C_{\text{crit}}, \\
L_z &= (2z-1) \sin\theta \sqrt{C - a^2 \cos^2\theta (1 - \mathcal{E}^2)},
\end{aligned} 
\quad (13) 
$$

which puts our $(x, y, z)$-integration region in $[0, 1] \times [0, 1] \times [0, 1]$.

The equivalent for a general inclination of the radius of the marginally bound circular orbit is found by setting $\mathcal{E}_{\min}(r, \theta) = 1$. The dark matter density vanishes for $r \le r_{\min}(\theta)$, since no bound orbit exists within this region.

Given a spin parameter $\tilde{a}$ and an inclination $\theta$, we can run the MC to calculate the current density for any $r > r_{\min}(\theta)$ with these simplifications. We use two standard numerical routines, VEGAS [60] and MISER [61], as implemented in the GNU Scientific Library [62]. When evaluating the density at a particular point, the numerical integrator performs $10^7$ function calls to reach a relative accuracy at the level of $\lesssim 1\%$ across the entire $r, \theta$ range.

Sample results are displayed in Figs. 7 and 8, where we use the same value of $f(\mathcal{E}, C, L_z) = f_0 = \frac{0.3 \text{ GeV/cm}^3}{(2\pi(100 \text{ km/s})^2)^{3/2}} = 5.1 \times 10^8 \text{ GeV/cm}^3$ as in SFW to allow a direct comparison of our results. The density plot in Fig. 9 provides a pictorial illustration of the density distribution.

Once more we find that the density increases with the spin parameter $\tilde{a}$, and the spike gets closer to the hole.

The decrease of the density as we get away from the equatorial plane is to be expected: for instance, only orbits with $L_z = 0$ can cross the axis. Unlike the spherical case, in which the components $L_x$ and $L_y$ of the angular momentum are also conserved, this restriction effectively reduces the available phase space and, consequently, the density. This is useful since the calculation for orbits that cross the axis is simpler and provides us with a lower bound on the density

<figure><img src="image_9.png" /><figcaption>FIG. 8. Density anisotropy for $\tilde{a} = 0.8$. The spike is shown at different angles with respect to the black hole rotation axis: on axis $\theta = 0$ (blue, dot-dashed), $\theta = \pi/3$ (purple, dashed) and equatorial $\theta = \pi/2$ (red, solid).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 9. Dark matter density in the $r-\theta$ plane for a spin parameter $\tilde{a} = 0.8$, and a constant initial distribution function. The axis of the black hole points vertically, and $r/Gm$ is plotted from 0 to 15. The density is axisymmetric about the spin axis.</figcaption></figure>

everywhere. Note that, as can be seen from comparing Figs. 7 and 8, although lower than on the plane, the density on the axes is still boosted for $\tilde{a} = 0.8$ compared to the case of Schwarzschild.

Indeed, using the last substitution in Eq. (13) and setting $\sin\theta = 0$ allows us to perform the integrals over $z$ and $C$ explicitly, obtaining:

$$ J_t = \frac{4\pi f_0}{\sqrt{\Delta}} \int_{\mathcal{E}_{\min}}^{1} d\mathcal{E} \sqrt{C_{\max}(\mathcal{E}) - C_{\text{crit}}(\mathcal{E})}, \quad (14) $$

where, as usual, $C_{\max}$ is obtained from $V(r) = 0$ and $C_{\text{crit}}$ is the critical Carter constant. Because $L_z = 0$, these are much
simpler constraints than in the general case: $C_{max}$ can be found algebraically and $C_{crit}$ requires numerically solving a straightforward nonlinear equation. On the axis, $J_{\phi} = 0$, and the density is easily found from Eq. (9) with $\Omega = 0$. This alternative route to calculate the density along the rotation axis of the black hole provides us with a useful check of the full MC evaluation.

## IV. GROWTH FROM A CUSPY HERNQUIST PROFILE

The constant distribution is an adequate proxy for a cored profile [20], but DM-only simulations tend to favor cuspy distributions, and we would like to know if the calculated effect is still sizeable for a cuspy profile. Following SFW, we consider a Hernquist profile [63]:

$$ \rho_H = \frac{\rho_0}{(r/r_s)(1+r/r_s)^3}, \quad (15) $$

which generates the Newtonian gravitational potential

$$ \Phi_H = -\frac{GM}{r_s + r}. \quad (16) $$

Here $\rho_0$ and $r_s$ are scale factors related to the total dark matter mass in the halo by $M = 2\pi\rho_0 r_s^3$, which we take to be $M = 10^{12} M_{\odot}$ for the Milky Way. The ergodic distribution function associated with the Hernquist profile can be found analytically [63]:

$$ f_H(\tilde{e}) = \frac{M}{\sqrt{2(2\pi)^3(GMr_s)^{3/2}}} \tilde{f}_H(\tilde{e}), \quad (17) $$

with

$$ \tilde{f}_H(\tilde{e}) = \frac{\sqrt{\tilde{e}}}{(1-\tilde{e})^2} \left[ (1-2\tilde{e})(8\tilde{e}^2 - 2\tilde{e} - 3) + \frac{3\sin^{-1}\sqrt{\tilde{e}}}{\sqrt{\tilde{e}(1-\tilde{e})}} \right], \quad (18) $$

and we have introduced a new dimensionless relative energy $\tilde{e} = -r_s E/GM$, which is related to the relativistic energy $\mathcal{E}$ per unit particle mass by

$$ \tilde{e} \equiv \frac{r_s}{GM} (1 - \mathcal{E}). \quad (19) $$

The halos found in simulations can be better fit by an NFW profile [2] with strongly correlated $\rho_0$, $r_s$. As a result, halos are essentially members of a one-parameter family. For a galactic mass halo, $r_s = 20$ kpc, and $\rho_0$ is then fixed by the total mass. We choose to work with a Hernquist profile with the same parameters because, as mentioned above, the distribution function can be calculated analytically, and

<figure><img src="image_11.png" /><figcaption>FIG. 10. Dark matter density in the equatorial plane with increasing Kerr parameter for an initial Hernquist distribution. The different lines show the Schwarzschild calculation from SFW (blue, dot-dashed), $\tilde{a} = 0.5$ (purple, dashed), and $\tilde{a} = 0.8$ (red, solid).</figcaption></figure>

both halos have the same cuspy $\propto 1/r$ behavior in the inner region $r \lesssim r_s$ giving rise to the same Newtonian spike [49].

Thus the only numerical difficulty introduced by this distribution is the evaluation of the radial action.

Using Eq. (16) in Eq. (10), we can write the radial invariant for the Hernquist potential as

$$ I_H^r = 2\sqrt{GMr_s} \int_{x_-}^{x_+} dx \left( \frac{2}{1+x} - 2\tilde{e} + \frac{\tilde{L}^2}{x^2} \right)^{1/2}, \quad (20) $$

where we introduced the dimensionless quantities $\tilde{L} = L\sqrt{GMr_s}$, $x = r/r_s$; and $x_{\pm}$ are the turning points of the orbit. As discussed in Sec. II B, $\tilde{L}$ and $\tilde{e}$ are obtained through the matching of Eqs. (10) to (11) for each point in the Kerr phase space ($\mathcal{E}$, $C$, $L_z$) with the correspondence Eq. (19).

As $\mathcal{E} \to 1$, the orbit becomes unbound and the corresponding $\tilde{e} \to 0$. In this limit, both radial invariants diverge as the upper turning point goes to infinity. To prevent numerical instabilities in the evaluation of the integrals, we follow SFW and remap the interval $[x_-, x_+]$ to [0, 1] at each step of the bisection method that is used to evaluate Eq. (20). We do a similar remapping to evaluate the Kerr invariant.<fnref n="3" />

The Hernquist distribution has most of its contribution to the density coming from the more deeply bound orbits, which are those that have $\tilde{e}$ closer to 1. Since we have argued in Sec. III that the presence of a deeper potential well for corotating orbits is what drives the enhancement of the spike, we expect it to increase more rapidly with increasing Kerr parameter. This is seen in Figs. 10 and 11 below: the peak of the spike with $\tilde{a} = 0.8$ is approximately 35% higher than the one with $\tilde{a} = 0.5$, whereas the corresponding boost for the constant

<footnote n="3">³Note that the Schwarzschild limit of the radial invariant in Eq. (3.19) in SFW is missing a factor of $r^2/\Delta = 1/(1 - 2Gm/r)$.</footnote>
<figure><img src="image_12.png" /><figcaption>FIG. 11. Density anisotropy for an initial Hernquist profile and $\tilde{a} = 0.8$. The spike is shown at different angles with respect to the black hole rotation axis: on axis $\theta = 0$ (blue, dot-dashed), $\theta = \pi/3$ (purple, dashed) and equatorial $\theta = \pi/2$ (red, solid).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 12. Dark matter density in the $r-\theta$ plane for a spin parameter $\tilde{a} = 0.8$, and an initial distribution function corresponding to a Hernquist profile.</figcaption></figure>

distribution is around 20%. Figure 12 shows the density distribution resulting from an initial Hernquist profile.

V. DISCUSSION

The spatial extension of the spike is of the order of the radius of gravitational influence of the black hole, $r_h = Gm/\sigma_v^2 \sim 1.7$ pc for the galactic center black hole, and the mass of the dark matter contained in this small volume is negligible compared to the mass of the black hole. Hence, as pointed out by SFW, the additional precession rate induced by the spike on the orbits of stars in the central cluster will be small, and this conclusion remains valid for a rotating black hole. As shown in Fig. 5 in SFW, the rapid fall off with distance of the effects of frame dragging

and departure from spherical symmetry makes their contribution to the pericenter advance of the orbits of stars subdominant when compared to the monopole Schwarzschild component over the range of semi-major axis $\gtrsim 0.1$ pc. The enhanced dark matter density due to the rotation of the hole, will neither interfere with a test of the black hole no-hair theorem using hypothetical stars with semimajor axes $\lesssim 0.2$ mpc [64], nor change the conclusions of a future experiment sensitive to precession rates at the level of $10 \ \mu\text{arcsec}$ per year that could discover the perturbing effects of the dark-matter for S2-type stars, with semimajor axes $\sim 10$ pc [21].

On the other hand, the enhanced density in the Kerr geometry could significantly alter the fluxes of high-energy radiation from dark matter annihilations in the central regions surrounding the black hole. A prediction of the effects of rotation depends on the underlying particle physics model and on the formation history of the super-massive black hole (see e.g. [65] for a review).

As mentioned above, we are assuming that the growth of the black hole is adiabatic, but several dynamical effects could weaken the density spike around the supermassive black hole. Although these effects are important, our main purpose is to understand the general relativistic effects close to the rotating black hole, and the following considerations provide an order of magnitude estimate in a sample of representative scenarios.

To estimate the effects of the spike on the dark matter annihilation fluxes we compare the line of sight integral for a beam of opening angle $\theta$ towards the Galactic Center:

$$J(\theta) = \frac{\langle \sigma v \rangle}{4\pi m_{\chi}^2} \int_0^{2\pi} d\phi \int_0^{\theta} d\psi \cos\psi \sin\psi \int ds \rho^2(r, \vartheta), \qquad (21)$$

with or without the presence of the spike. In Eq. (21), $r = \sqrt{R^2 + s^2 - 2Rs \cos\psi}$ is the Boyer-Lindquist coordinate, $s$ is the radial coordinate from the Earth to the annihilation point, $R = 8.5$ kpc is the distance from the Earth to the Galactic Center, and $\vartheta$ is the angle relative to the equatorial plane of the black hole.

For a thermal relic, $\langle \sigma v \rangle = 3 \times 10^{-26} \text{ cm}^3/\text{s}$, of mass $m_{\chi} = 100$ GeV, we find a flux of $6.3 \times 10^{-9} \text{ cm}^{-2} \text{s}^{-1}$ for an opening angle of $1^\circ$ from annihilations in the halo with only the underlying Hernquist profile included.

To evaluate the integral in Eq. (21), we fit the profiles that were previously calculated using MC techniques to profiles of the form

$$\rho(r, \vartheta) = \frac{A}{x^p} \left( 1 - \frac{r_{\min}(a, \vartheta)}{r} \right)^n, \qquad (22)$$

with all coefficients being allowed to vary with $\vartheta$. The expression in Eq. (22) is a generalization of Eq. (9) in [20]
<table><caption>TABLE I. Boost factors for different Kerr spin parameters $\tilde{a}$ for the full spike with no annihilation normalized to the Schwarzschild spike.</caption><thead><tr><th>$\tilde{a}$</th><th>$J(a)/J(0)$</th></tr></thead><tbody><tr><td>0</td><td>1</td></tr><tr><td>0.5</td><td>1.11</td></tr><tr><td>0.6</td><td>1.14</td></tr><tr><td>0.7</td><td>1.22</td></tr><tr><td>0.8</td><td>1.38</td></tr><tr><td>0.9</td><td>1.59</td></tr><tr><td>0.998</td><td>1.97</td></tr></tbody></table>

to allow for a $\theta$-dependent end point for the spike. This expression is then matched to the power law $B/x^{\gamma_{sp}}$ using smooth functions to improve the fit and to give a reasonable estimate of $\rho$. We use $\gamma_{sp} = 2.33$ corresponding to the Newtonian spike generated by a $1/r$ NFW or Hernquist cusp.

We extend our spike profile until the density is equal to that of the underlying Hernquist profile. This happens at 12.4 pc, and we take the fiducial $1/r$ halo shape beyond this point. Putting all the pieces together, our model for the central part of the halo is equivalent to the canonical model used in [14], aside from the fact that we are assuming an underlying Hernquist profile with $r_s = 20$ kpc and the extension of the spike. For the Schwarzschild geometry, the spike enhances the flux by a factor of $1.93 \times 10^9$ relative to the initial Hernquist profile. In Table I, we present the flux enhancement normalized to the Schwarzschild one, as that is less sensitive to the normalization of the underlying density profile.

Note that, as mentioned above, the influence radius of the black hole is $r_h \sim 1.7$ pc, so at these distances the spike should start being modified. Also, our calculation of the Kerr spike does not consistently take into account the gravitational field generated by the dark matter distribution itself. We leave for future work a proper treatment of this effect along the lines of the grid calculation in [49] for the nonrelativistic case. Let us stress that it is only in the transition region between the black hole dominated field and the smooth underlying halo that we expect significant changes. Moreover, the relative boost factors quoted in Table I will remain mostly unchanged.

Our expression for the spike does not include the effects of dark matter annihilations which will deplete and weaken the density profile. SFW followed the strategy used in [20], and considered a constant core within the radius $r_{ann}$ determined by the location where the density equals $\rho_{ann} = m_{\chi}/\sigma v t_{bh}$, which for our thermal relic turns out to be $r_{ann} = 1.3 \times 10^{-2}$ pc, assuming that the annihilation process has been acting over $t_{bh} = 10^{10}$ yr. It was pointed out in [42] that a plateau is in equilibrium if all the dark matter particles move in circular orbits, and the density forms a “weak cusp” $\propto r^{-1/2}$ for a more realistic isotropic distribution of DM velocities. This behavior has been

confirmed by integrating the Boltzmann equation in [43] for the case of s-wave annihilation. If we assume that the DM is distributed as $r^{-1/2}$ within $r_{ann}$ then the density is sufficiently low that the spike goes away. The effects of the black hole spin are therefore washed out and we get the same boost factor for a Schwarzschild as for a nearly extreme $\tilde{a} = 0.998$ black hole. In both cases the profile close to the hole is the same $r^{-1/2}$ weak spike, and the ratio of the flux coming from the spike to that coming from the smooth halo is $2.8 \times 10^3$. This large enhancement is due to the spatial extension of our spike, which shows that most of the annihilation signal comes from the outer regions of the spike, where relativistic corrections are not important.

The weak $r^{-1/2}$ is not generic, but depends on the s-wave nature of the annihilating process. For p-wave annihilation, the final cusp is even weaker and again erases any effects due to the Kerr spin. On the other hand, a possible detection of an identified x-ray line at $E \approx 3.55$ keV from the Perseus cluster and the Andromeda galaxy [66,67] have been attributed to a sterile neutrino [68]. This dark matter candidate decays emitting an x-ray line and its density profile is the Kerr spike without any attenuation due to self-annihilation. The x-ray flux due to dark matter decays, however, is proportional to the density as opposed to the density squared behavior of the annihilation signal shown in Eq. (21). Much like in the case of the precession rates induced by the spike on the orbits of stars, we expect the additional decay flux induced by the spike to be small, and the effects due to the rotation of the black hole to be subdominant.

Self-interacting dark matter [69] with a cross section as large as $\sigma/m_{\chi} \approx 0.1 \text{ cm}^2/\text{g}$ has been invoked to address several discrepancies between numerical predictions of cold dark matter models and observations of subgalactic scale structures [70,71]. The effects of self-interactions have been shown in [72] to replenish the weak cusp, giving rise to a steeper profile depending on the velocity dependence of the cross section. Since a detailed study of the final relaxed distribution using the exact phase-space distribution of the Kerr spike is missing, we display in Table I the radiation fluxes as arising from a full Kerr spike with no weakening due to annihilation. The results show that when going from a nonrotating to an almost extremal black hole, the boost factor is almost doubled.

## VI. CONCLUSIONS

We have extended the analysis performed in SFW to include the effects of black hole spin. Our findings show that the spike persists around a rotating black hole and, furthermore, that it is enhanced.

Since the total mass contained in the spike is not very large, effects that depend on the total mass of the spike such as the stellar precession studied in [21] or fluxes from decaying dark matter will essentially remain unaltered by the inclusion of rotation. Our results are summarized in
<figure><img src="image_14.png" /><figcaption>FIG. 13. Equatorial (solid blue) and on-axis (solid red) density obtained for the near-extreme black hole, compared with previous calculations assuming a Schwarzschild hole (dot-dashed purple), and the nonrelativistic estimate in [20] (dashed green). The initial DM distribution before the growth of the black hole is a Hernquist profile shown as the dotted red line.</figcaption></figure>

Fig. 13, which shows the largest density enhancement, obtained for a near-extreme black hole (ā = 0.998) in the equatorial plane, as well as the density along the spin axis, comparing them to previous calculations.

However, the large growth of the spike could have consequences for observables related to dark matter annihilation, which depend on the density squared. Further work remains to be done here to properly implement the evolution of the spike and the presence of a weak cusp [42,43], which are both necessary for the extraction of a gamma-ray signal.

Although we have focused our attention on the signals from the Galactic Center, black holes are ubiquitous in nature. Our findings could potentially impact the cumulative effects of dark matter spikes on the diffuse gamma-ray background [47], or the signals from local dwarf spheroidal galaxies [73].

## ACKNOWLEDGMENTS

We are grateful to Tongbo Liu, Oleg Ruchayskiy, Laleh Sadeghian, and Jessie Shelton for useful discussions. F. F. acknowledges the hospitality of the Aspen Center for Physics, which is supported by the National Science Foundation. F. F. and A. M. R. were supported in part by the U.S. Department of Energy, Office of High Energy Physics, under Award No. DE-FG02-91ER40628. C. M. W. was supported in part by the National Science Foundation under Grant No. 16-00188.

## APPENDIX A: ADIABATIC INVARIANCE OF THE DISTRIBUTION FUNCTION

We here prove the important property that the distribution function is an adiabatic invariant, which we use in our code. A proof for nonrelativistic spherical systems is given in Young [48], which was generalized to the general

relativistic Schwarzschild case by Sadeghian [58]. The extension to the Kerr geometry uses the same ideas, but the change of variables that must be performed in the adiabatic evolution is two-dimensional instead of one-dimensional.

We start by integrating the current density Eq. (8) to find the mass enclosed in a hypersurface of constant time. The future-pointing normal vector to the surface is $n_α = -(-g^{tt})^{-1/2}∂_α t$, and the three-dimensional surface element can be written as [74]:

$$d^3S_α = -\delta_α^0 \frac{\sqrt{g_S}}{\sqrt{-g^{tt}}} d^3x, \quad (A1)$$

where

$$g_S = \frac{\Sigma^4}{\Delta} g_{\phi\phi}, \quad (A2)$$

is the determinant of the metric induced on the hypersurface.

The enclosed mass is therefore:

$$\begin{aligned} M &= -\int d^3S_α J^α \\ &= \int dr dθ dφ \frac{\sqrt{g_S}}{\sqrt{-g^{tt}}} (g^{tt} J_0 + g^{tφ} J_φ) \\ &= 2 \int \frac{dr dθ dφ dE dC dL_z}{r^2 \sin θ \sqrt{V(r)} \sqrt{U(θ)} \sqrt{-\Delta g^{tt}}} \frac{\Sigma^2 \sqrt{g_{φφ}}}{\sqrt{-\Delta g^{tt}}} \\ & \quad \times (g^{tφ} L_z - g^{tt} E) f(E, C, L_z). \end{aligned} \quad (A3)$$

Interchanging the order of integration, the integrals over coordinates give us a distribution $N(E, C, L_z)$ of particles per unit conserved quantity. Using the relation $g^{tt} = -g_{\phi\phi}/(\Delta\sin^2\theta)$ for Kerr, we obtain:

$$N(E, C, L_z) = 4\pi \int \frac{dr d\theta}{r^2 \sqrt{U(\theta)} \sqrt{V(r)}} \Sigma^2 \\ \times (g^{t\phi} L_z - g^{tt} E) f(E, C, L_z). \quad (A4)$$

In the above expression, the limits of integration are the radial turning points of $V(r)$ and the angular turning points of $U(θ)$. Under adiabatic evolution, the constants of the motion will change to new values $E^*, C^*, L_z^*$, in such a way that $N(E, C, L_z)dE dC dL_z = N^*(E^*, C^*, L_z^*)dE^* dC^* dL_z^*$. The new values of the constants of the motion are determined by the invariance of the action integrals in Eq. (11). If we assume that the initial and final states are axisymmetric, then the invariance of $I_φ$ gives $L_z = L_z^*$. The remaining two integrals allow us to relate $N$ and $N^*$ through the chain rule of multivariable calculus,
$$ \frac{\partial(\mathcal{E}^*, C^*)}{\partial(\mathcal{E}, C)} = \frac{\partial(\mathcal{E}^*, C^*)}{\partial(I_r^*, I_\theta^*)} \frac{\partial(I_r, I_\theta)}{\partial(\mathcal{E}, C)}. $$

Taking partial derivatives in Eq. (11), we find:

$$ \begin{aligned} \mathfrak{F} &\equiv \frac{\partial(I_r, I_\theta)}{\partial(\mathcal{E}, C)} \\ &= 2 \int \frac{\mathrm{d}r \mathrm{d}\theta \Sigma^2(\mathcal{E}g_{\phi\phi} + L_z g_{t\phi})}{r^2 \sin^2\theta \Delta \sqrt{U(\theta)} \sqrt{V(r)}}. \end{aligned} \quad (A5) $$

Using $g^{tt} = -g_{\phi\phi}/(\Delta\sin^2\theta)$, $g^{t\phi} = g_{t\phi}/(\Delta\sin^2\theta)$ and comparing to Eq. (A4), we immediately obtain

$$ N(\mathcal{E}, C, L_z) = 2\pi \mathfrak{F} f(\mathcal{E}, C, L_z). $$

Therefore, we find

$$ \begin{aligned} N(\mathcal{E}, C, L_z) \mathrm{d}\mathcal{E} \mathrm{d}C \mathrm{d}L_z \\ &= 2\pi \mathfrak{F} f(\mathcal{E}, C, L_z) \mathrm{d}\mathcal{E} \mathrm{d}C \mathrm{d}L_z \\ &= 2\pi \mathfrak{F}^* f^*(\mathcal{E}^*, C^*, L_z^*) \mathrm{d}\mathcal{E}^* \mathrm{d}C^* \mathrm{d}L_z^* \\ &= 2\pi \mathfrak{F}^* f^*(\mathcal{E}^*, C^*, L_z^*) \\ & \quad \times \frac{\partial(\mathcal{E}^*, C^*)}{\partial(I_r^*, I_\theta^*)} \frac{\partial(I_r, I_\theta)}{\partial(\mathcal{E}, C)} \mathrm{d}\mathcal{E} \mathrm{d}C \mathrm{d}L_z \\ &= 2\pi \mathfrak{F}^* f^*(\mathcal{E}^*, C^*, L_z^*) \frac{\mathfrak{F}}{\mathfrak{F}^*} \mathrm{d}\mathcal{E} \mathrm{d}C \mathrm{d}L_z \\ &= 2\pi \mathfrak{F} f^*(\mathcal{E}^*, C^*, L_z^*) \mathrm{d}\mathcal{E} \mathrm{d}C \mathrm{d}L_z, \end{aligned} $$

which demonstrates the adiabatic invariance of the distribution function,

$$ f(\mathcal{E}, C, L_z) = f^*(\mathcal{E}^*, C^*, L_z^*). $$

## APPENDIX B: PHASE SPACE FOR ORBITS THAT REMAIN ON THE EQUATORIAL PLANE

We will set $G = 1$ and work with the dimensionless quantities:

$$ \tilde{a} \equiv \frac{a}{m} \qquad \tilde{L}_z \equiv \frac{L_z}{m} \qquad \tilde{C} \equiv \frac{C}{m^2} \qquad x \equiv \frac{r}{m}, $$

and we will drop the tildes for the rest of this appendix. The black hole horizon is located at

$$ x_{\text{horizon}} = 1 + \sqrt{1 - a^2}, \quad (B1) $$

and the boundary of the ergosphere is

$$ x_{\text{ergosphere}} = 1 + \sqrt{1 - a^2 \cos^2 \theta}, \quad (B2) $$

which is equal to $x_{\text{ergosphere}} = 2$ for the equatorial $\theta = \pi/2$ latitude.

For planar trajectories $L_z = \pm\sqrt{C}$, with the plus (minus) sign corresponding to corotating (counter-rotating) orbits. Hence, the subset of orbits that we are considering is defined by two quantities $\mathcal{E}$ and $L_z$. These are constrained to satisfy $\mathcal{E} \le 1$, which bounds the energy from above; $V(r) \ge 0$, which results in an upper bound on the angular momentum $L_z^{\text{max}}$; and the capture condition

$$ V = 0 = \frac{\mathrm{d}V}{\mathrm{d}x}, \quad (B3) $$

which will determine the critical angular momentum, $L_z^{\text{crit}}$, and the minimum energy, $\mathcal{E}_{\text{min}}$. Note that the constraint $U(\theta) \ge 0$ is trivially satisfied.

The effective potential for orbits in the equatorial plane reads:

$$ V|_{C=L_z^2} = \frac{\mathcal{E}^2(x^3 + a^2(2+x)) - 4aL_z\mathcal{E} + (2-x)(L_z^2 + x^2) - a^2x}{x^3}. \quad (B4) $$

We will use this expression to obtain $L_z^{\text{max}}$, $L_z^{\text{crit}}$ and $\mathcal{E}_{\text{min}}$, at a fixed distance $x$ from the hole.

### 1. $L_z^{\text{max}}$

For a fixed energy $\mathcal{E} < 1$, the constraint $V \ge 0$ puts an upper bound on $L_z$. We can explicitly find the boundary $V = 0$, since $V$ is quadratic in $\mathcal{E}$ and $L_z$:

$$ L_z^{\text{max}} = \frac{-2a\mathcal{E} \pm \sqrt{x(a^2 - x(\mathcal{E} - 2))(2 - (1 - \mathcal{E}^2)x)}}{x - 2}. \quad (B5) $$

We generally obtain two solutions: a positive value of $L_z$, which is the maximum for corotating orbits, and a negative

one corresponding to the smallest $L_z$ for counter-rotating orbits (which we also name $L_z^{\text{max}}$). But these solutions might not exist for all $x$, since the polynomial inside the square root can become negative. Being a quartic polynomial in $x$ with negative quartic coefficient, $\mathcal{E}^2 - 1$, it will become negative (i.e. no $L_z^{\text{max}}$) at large $x$. This reflects the fact that far away orbits are closer to being unbound as discussed below [cf. Eq. (B8)].

More specifically, the discriminant in Eq. (B5) vanishes at the nonphysical points $x = 0$ and $x = 1 - \sqrt{1 - a^2}$, which are inside the horizon. It also vanishes at

$$ x = x_{\text{horizon}} \quad \text{and} \quad x = \frac{2}{1 - \mathcal{E}^2} > 2, \quad (B6) $$
so the square root is well defined between these two locations and we will find two values of $L_z^{\text{max}}$.

This is true for points outside of the ergosphere, $x > 2$. At the ergosphere, Eq. (B4) is linear in $L_z$, and there is only one solution. A closer look shows that the positive root smoothly tends to this single root,

$$L_z^{\text{max}} = \frac{2a^2\mathcal{E}^2 + 4\mathcal{E}^2 - a^2}{2a\mathcal{E}},$$

which may be positive.

The negative branch in Eq. (B5) diverges as $x \to 2$. This does not mean that counter-rotating orbits with arbitrarily large $|L_z|$ are allowed, since we will see that they cease to exist that close to the hole. For the Schwarzschild case SFW found that the critical value of the angular momentum limited orbits to lie beyond the unstable marginally bound orbit at $x = 4$ for $\mathcal{E} = 1$. The minimum energy $\mathcal{E} = \sqrt{8/9}$ is attained at the location of the ISCO, $x = 6$. For a Kerr hole the locations of these orbits get closer to the horizon for corotating orbits ($x = 1$ for the extremal $a = 1$ case) and pushed away for counter-rotating ones ($x = 5.83$ and $x = 9$). As we prove below, there are no orbits within the radius of the marginally bound orbit,

$$x_{\text{mb}} = 2 \mp a + 2\sqrt{1 \mp a}, \qquad (\text{B7})$$

and we will not have to worry about regions where Eq. (B5) does not provide a valid bound. Nevertheless, there might be corotating bound orbits within the ergosphere since Eq. (B7) can be $< 2$.

$$ \begin{aligned} V(x, \mathcal{E}, L_z^{\text{crit}}) = 0 &\Rightarrow \\ 0 &= -9(L_z^{\text{crit}} - a\mathcal{E})^2 \\ &\times [(\mathcal{E}^2 - 1)L_z^{\text{crit}6} + (36\mathcal{E}^2 - 27\mathcal{E}^4 - 3a^2(1 - \mathcal{E}^2)^2 - 8)L_z^{\text{crit}4} + 36a\mathcal{E}(2 - 5\mathcal{E}^2 + 3\mathcal{E}^4)L_z^{\text{crit}3} \\ &+ (2a^2(10 - 91\mathcal{E}^2 + 162\mathcal{E}^4 - 81\mathcal{E}^6) + 3a^4(\mathcal{E}^2 - 1)^3 - 16)L_z^{\text{crit}2} + 4a\mathcal{E}(8 + 9a^2(1 - \mathcal{E}^2)^2(-1 + 3\mathcal{E}^2))L_z^{\text{crit}1} \\ &- 16a^2\mathcal{E}^2 - a^6(1 - \mathcal{E}^2)^4 - a^4(1 - \mathcal{E}^2)^2(-1 - 18\mathcal{E}^2 + 27\mathcal{E}^4)], \end{aligned} \qquad (\text{B10}) $$

which contains a sixth-order polynomial in $L_z^{\text{crit}}$ that can be numerically solved for given $a$ and $\mathcal{E}$. Since $L_z^{\text{crit}} \le L_z^{\text{max}}$, we select the largest of all the roots that are smaller than $L_z^{\text{max}}$.

The value of $L_z^{\text{crit}}$ will depend on our choice of $\mathcal{E}$. For corotating orbits it can be shown from Eq. (B5) that $L_z^{\text{max}}$ is a monotonically increasing function of $\mathcal{E}$, so it will be largest for $\mathcal{E} = 1$, while $L_z^{\text{crit}}$ is a monotonically decreasing function of $\mathcal{E}$. On the other hand, keeping the value of $L_z$

Note also that, for a given $x$, there is a minimum value of the energy since, for $V$ to be positive, we need $x$ to be to the left of the second root in Eq. (B6), so

$$\mathcal{E}^2 \ge 1 - \frac{2}{x}. \qquad (\text{B8})$$

## 2. $L_z^{\text{crit}}$ and $\mathcal{E}_{\text{min}}$

Let us examine the capture condition, which is going to provide us with $L_z^{\text{crit}}$ and $\mathcal{E}_{\text{min}}$. The first equality in Eq. (B3) is a cubic equation in $x$, while the second requires solving a 2nd order polynomial.

The solution to the second equation is

$$ \begin{aligned} \frac{dV}{dx} = 0 &\Rightarrow x_{\mp} = (a^2(1 - \mathcal{E}^2) \\ &+ L_z^2 \mp \sqrt{(L_z^2 + a^2(1 - \mathcal{E}^2))^2 - 12(L_z - a\mathcal{E})^2})/2. \end{aligned} \qquad (\text{B9}) $$

The smaller root $x_{\mp}$ will correspond to the turning point when the constraint $V = 0$ is also satisfied. It has the correct Schwarzschild limit:

$$\lim_{x \to a} x_{\mp} = \frac{L_z^2 - L_z \sqrt{L_z^2 - 12}}{2} = \frac{6}{1 + \sqrt{1 - 12/L_z^2}},$$

which is the location of the unstable circular orbit.

We can now plug this value back in the equation $V = 0$ and solve for $L_z^{\text{crit}}$:

constant, there is a minimum energy $\mathcal{E}_{\text{min}}$, which satisfies $V = 0$.

This allows us to find the boundary of momentum-space at the location $x$ using the following numerical scheme depicted in Fig. 14. Starting with $\mathcal{E} = 1$,

(i) Find $L_z^{\text{max}}$ from Eq. (B5), and $L_z^{\text{crit}}$ solving Eq. (B10). As long as $L_z^{\text{max}} \ge L_z^{\text{crit}}$, the volume in momentum space at $x$ is nonzero.

(ii) Solve $V(L_z = L_z^{\text{crit}}) = 0$ to find $\mathcal{E}_{\text{min}}$.
<figure><img src="image_15.png" /><figcaption>FIG. 14. Momentum-space at the location $x = 10$, which is shown as a red dot. Top left: Starting from allowed values of $\mathcal{E} = \mathcal{E}_0 = 0.9997$ and $L_z = 4$, the red dot is well within the allowed region. Top right: For this value of $\mathcal{E}$, changing $L_z \to L_z^{\text{max}} = 4.829$ from Eq. (B5) pulls the effective potential $V_{\text{eff}} \equiv -V(r)$ up and brings the red dot to the edge of the allowed region. Bottom left: Keeping $\mathcal{E} = \mathcal{E}_0$ but changing $L_z \to L_z^{\text{crit}} = 2.9$ brings the effective potential down, so that the orbit is almost captured. Bottom right: Once we have the $L_z$ bounds, lowering the energy while keeping $L_z = L_z^{\text{crit}}$ will pull the effective potential curve up. When it touches the red dot we have found $\mathcal{E}_{\text{min}} = 0.938$. Setting $\mathcal{E}_0 = \mathcal{E}_{\text{min}}$ the process can be iterated to find $L_z^{\text{max}} = 2.99$, $L_z^{\text{crit}} = 2.73$, and the next $\mathcal{E}_{\text{min}} = 0.931$.</figcaption></figure>

Note that we have found the minimum energy of all the orbits with $L_z = L_z^{\text{crit}}$, where $L_z^{\text{crit}}$ was found by setting $\mathcal{E} = 1$. But we can now set $\mathcal{E} = \mathcal{E}_{\text{min}}$ and repeat the steps above to find a smaller $L_z^{\text{max}}$ and a larger $L_z^{\text{crit}}$; and a smaller $\mathcal{E}_{\text{min}}$. This process can be iterated until we get to an energy for which $L_z^{\text{max}} = L_z^{\text{crit}}$, this is the minimum energy for any $L_z$ allowed at the location $x$. An analogous strategy can be followed for counter-rotating orbits.

As we get closer to the hole, the phase space gets reduced. Orbits with smaller values of $\mathcal{E}$, which were allowed at larger distances, are now pulled in and trapped. Eventually, we reach a point where only a single orbit with $\mathcal{E} = 1$ is allowed. In this case, the iteration process above will end immediately, since $L_z^{\text{max}} = L_z^{\text{crit}}$ for $\mathcal{E} = 1$ at that point. We can find this location by solving for $x$ from the equation $L_z^{\text{max}} = L_z^{\text{crit}}$ with $\mathcal{E} = 1$.

It turns out that this can be done analytically, since for $\mathcal{E} = 1$ Eq. (B10) simplifies to

$$0 = -9(a - L_z)^2(-16a^2 + 32aL_z - 16L_z^2 + L_z^4).$$

The roots of this polynomial can be found in closed form:

$$L_z^{\text{crit}} = a, 2(1 - \sqrt{1 - a}), \quad 2(1 + \sqrt{1 - a}), \\ 2(-1 - \sqrt{1 + a}), \quad 2(-1 + \sqrt{1 + a}). \qquad (\text{B11})$$

On the other hand, Eq. (B5) becomes for $\mathcal{E} = 1$:

$$L_z^{\text{max}} = \frac{-2a \pm \sqrt{2\sqrt{x(a^2 + (x-2)x)}}}{x-2}. \qquad (\text{B12})$$

We can solve for the values $x$ where $L_z^{\text{crit}} = L_z^{\text{max}}$ by equating each of the roots in Eq. (B11) to each of the two branches in Eq. (B12). Equating the negative branch of Eq. (B12) to the only negative root in Eq. (B11) and solving for $x$ we find for counter-rotating orbits,

$$x_{\text{min}}^{L_z<0} = 2 + a + 2\sqrt{1 + a}.$$

For corotating orbits, we equate the four positive roots in Eq. (B11) to the positive branch of Eq. (B12). The solution for $L_z^{\text{max}} = a$ occurs at $x = 0, a^2/2$, which are inside the horizon and thus nonphysical; $L_z^{\text{max}} = 2(\pm 1 \mp \sqrt{1 \mp a})$ are located at $x = 2 \pm a - 2\sqrt{1 \pm a}$, which are also inside the horizon. The only physical solution is
<figure><img src="image_16.png" /><figcaption>FIG. 15. Comparison of the edges of the phase-space volume as a function of distance $x = r/m$ from the hole for corotating orbits (solid blue), counter-rotating orbits (dashed magenta), and Schwarzschild (dot-dashed red). The left panel shows $\mathcal{E}_{\min}$, while the range of allowed values of $L_z$ is shown in the right one as the regions bounded by $L_z^{\max}$ and $L_z^{\text{crit}}$. The Kerr hole has spin $a = 0.75$, and the sign of the angular momentum for the counter-rotating orbits has been flipped. The dotted vertical lines mark the location of the marginally bound orbit, where the phase-space vanishes.</figcaption></figure>

$$x_{\min}^{L_z>0} = 2 - a + 2\sqrt{1-a}.$$

The locations above coincide with Eq. (B7). Hence, the volume of momentum-space available vanishes at the location of the marginally bound orbit, and only corotating orbits exist inside the ergosphere.

The boundaries derived above are shown in Fig. 15, which also shows the Schwarzschild limit for comparison. Note that in the latter case, the absolute minimum energy $\mathcal{E}_{\min} = \sqrt{8/9} = 0.94$ is attained at $x = 6$ in agreement with the results in SFW.

### 3. Density on the equatorial plane

Let us consider the two-dimensional equatorial plane described by the restriction of the Kerr geometry to $\theta = \pi/2$. The definition of $J^\mu$ in Eq. (1) is still valid, but now the indices run over $\mu = t, r, \phi$ only. Our geodesics are described by the constants of motion $\mathcal{E}, L_z, \mu$, and the Jacobian of the change is

$$\mathcal{J} = \left| \frac{\partial(\mathcal{E}, L_z, \mu)}{\partial(p^t, p^r, p^\phi)} \right| = \frac{\Delta|u_r|\sin^2\theta}{\mu^2}, \quad (B13)$$

where we are temporarily using dimensionful quantities to parallel the discussion in the main text. Now there is only an additional factor of 2 to account for the sign in $p^\phi$, and $\sqrt{-g} = r$. Putting things together and using $|u_r| = r^2/\Delta\sqrt{V(r)}$, we obtain

$$\sqrt{-g} d p^t d p^r d p^\phi = \frac{2\mu^2}{r\sqrt{V(r)}} d\mathcal{E} dL_z d\mu.$$

We have the relation $f^{(3)}(x, p) = \mu^{-2} f(\mathcal{E}, L_z)\delta(\mu - \mu_0)$, where $f(\mathcal{E}, L_z)$ has units of surface mass density. Considering a constant distribution on the plane $f(\mathcal{E}, L_z) = f_{\text{eq}}$, we obtain

$$J_{\mu} = \int f^{(3)}(x, p) u_{\mu} \sqrt{-g} d^3 p = \frac{2 f_{\text{eq}}}{r} \int u_{\mu} \frac{d\mathcal{E} dL_z}{\sqrt{V(r)}}. \quad (B14)$$

Going back to using dimensionless quantities, e.g. $L_z = Gm\tilde{L}_z$ and dropping the tildes again, we obtain more explicitly

$$\begin{align}
J_t &= -2 f_{\text{eq}} \sqrt{x} \int_{\mathcal{E}_{\min}}^{1} \mathcal{E} d\mathcal{E} \mathfrak{F}_t(\mathcal{E}), \\
J_{\phi} &= 2 f_{\text{eq}} \sqrt{x} (Gm) \int_{\mathcal{E}_{\min}}^{1} d\mathcal{E} \mathfrak{F}_{\phi}(\mathcal{E}), \tag{B15}
\end{align}$$

where

$$\begin{align}
\mathfrak{F}_t(\mathcal{E}) &\equiv \int_{\Delta L_z} \frac{dL_z}{\sqrt{\tilde{V}(x)}} \\
\mathfrak{F}_\phi(\mathcal{E}) &\equiv \int_{\Delta L_z} \frac{L_z dL_z}{\sqrt{\tilde{V}(x)}}, \tag{B16}
\end{align}$$

and $\tilde{V}(x)$ is the numerator in Eq. (B4). The $\mathcal{E}$-dependent region of integration $\Delta L_z = (L_z^-, L_z^{\text{crit},-}) \cup (L_z^{\text{crit},+}, L_z^+)$ includes both corotating and counter-rotating orbits, and $L_z^+, L_z^-$ are the respective values of $L_z^{\max}$ in Eq. (B5).

The integral over $L_z$ can be found analytically by means of an Euler substitution [75]. The end result is
$$ 
\begin{aligned}
\mathfrak{F}_t(\mathcal{E}) &= \frac{2}{\sqrt{|x-2|}} \left[ \theta(x - x_{\text{mb}}^+) \left( \theta(x-2) \arctan \frac{1}{\sqrt{\kappa^+}} + \theta(2-x) \arctan \frac{1}{\sqrt{-\kappa^+}} \right) + \theta(x - x_{\text{mb}}^-) \arctan \frac{1}{\sqrt{\kappa^- - 1}} \right], \\
\mathfrak{F}_\phi(\mathcal{E}) &= \frac{1}{\sqrt{|x-2|}} \left[ \theta(x - x_{\text{mb}}^+) \left( \theta(x-2) \left\{ (L_z^+ + L_z^-) \arctan \frac{1}{\sqrt{\kappa^+}} + \sqrt{(L_z^+ - L_z^{\text{crit,+})}(L_z^{\text{crit,+}} - L_z^-)} \right\} \right. \right. \\
&\qquad \left. + \theta(2-x) \left\{ (L_z^+ + L_z^-) \arctan \frac{1}{\sqrt{-\kappa^+}} - \sqrt{(L_z^{\text{crit,+}} - L_z^+)(L_z^{\text{crit,+}} - L_z^-)} \right\} \right) \\
&\qquad \left. + \theta(x - x_{\text{mb}}^-) \left\{ (L_z^+ + L_z^-) \arctan \frac{1}{\sqrt{\kappa^- - 1}} - \sqrt{(L_z^+ - L_z^{\text{crit,-})}(L_z^{\text{crit,-}} - L_z^-)} \right\} \right],
\end{aligned} 
\quad (B17) $$

where we have defined

$$ \kappa^+ \equiv \frac{L_z^{\text{crit,+}} - L_z^-}{L_z^+ - L_z^{\text{crit,+}}} \qquad \kappa^- \equiv \frac{L_z^+ - L_z^-}{L_z^{\text{crit,-}} - L_z^-}. \quad (B18) $$

Since we are working with a geometry that has only two spatial dimensions, our final result is a surface mass density. Working instead with the full Kerr geometry and the restriction $f(\mathcal{E}, C, L_z) = f_{\text{eq}}\delta(u_\theta)$ on the distribution

function, results in a volume mass density, which is what we used in Sec. III. The only difference between the results of this appendix and those in Sec. III is a factor of $2/r$ between the currents in Eq. (B15) and those in Eq. (12). When calculating quantities integrated over space, such as masses or annihilation fluxes, the difference is only a factor of 2, since the three-dimensional Jacobian has an extra factor of $r$ compared to the planar geometry.

[1] F. Iocco, M. Pato, and G. Bertone, Nat. Phys. 11, 245 (2015).

[2] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996).

[3] G. R. Blumenthal, S. M. Faber, R. Flores, and J. R. Primack, Astrophys. J. 301, 27 (1986).

[4] O. Y. Gnedin, D. Ceverino, N. Y. Gnedin, A. A. Klypin, A. V. Kravtsov, R. Levine, D. Nagai, and G. Yepes, arXiv:1108.5736.

[5] S. Peirani, S. Kay, and J. Silk, Astron. Astrophys. 479, 123 (2008).

[6] S. Mashchenko, H. M. P. Couchman, and J. Wadsley, Nature (London) 442, 539 (2006).

[7] F. Governato, A. Zolotov, A. Pontzen, C. Christensen, S. H. Oh, A. M. Brooks, T. Quinn, S. Shen, and J. Wadsley, Mon. Not. R. Astron. Soc. 422, 1231 (2012).

[8] A. Di Cintio, C. B. Brook, A. V. Macci, G. S. Stinson, A. Knebe, A. A. Dutton, and J. Wadsley, Mon. Not. R. Astron. Soc. 437, 415 (2014).

[9] M. Kuhlen, M. Vogelsberger, and R. Angulo, Phys. Dark Universe 1, 50 (2012).

[10] D. Hooper, Phys. Dark Universe 15, 53 (2017).

[11] M. Benito, N. Bernal, N. Bozorgnia, F. Calore, and F. Iocco, J. Cosmol. Astropart. Phys. 02 (2017) 007.

[12] L. Goodenough and D. Hooper, arXiv:0910.2998.

[13] D. Hooper and T. Linden, Phys. Rev. D 84, 123005 (2011).

[14] B. D. Fields, S. L. Shapiro, and J. Shelton, Phys. Rev. Lett. 113, 151302 (2014).

[15] F. Calore, I. Cholis, C. McCabe, and C. Weniger, Phys. Rev. D 91, 063003 (2015).

[16] F. Calore, N. Bozorgnia, M. Lovell, G. Bertone, M. Schaller, C. S. Frenk, R. A. Crain, J. Schaye, T. Theuns, and J. W. Trayford, J. Cosmol. Astropart. Phys. 12 (2015) 053.

[17] M. Ackermann et al. (Fermi-LAT), Astrophys. J. 840, 43 (2017).

[18] C. Karwin, S. Murgia, T. M. P. Tait, T. A. Porter, and P. Tanedo, Phys. Rev. D 95, 103005 (2017).

[19] R. Genzel, F. Eisenhauer, and S. Gillessen, Rev. Mod. Phys. 82, 3121 (2010).

[20] P. Gondolo and J. Silk, Phys. Rev. Lett. 83, 1719 (1999).

[21] L. Sadeghian, F. Ferrer, and C. M. Will, Phys. Rev. D 88, 063522 (2013).

[22] M. Volonteri, P. Madau, E. Quataert, and M. J. Rees, Astrophys. J. 620, 69 (2005).

[23] E. Berti and M. Volonteri, Astrophys. J. 684, 822 (2008).

[24] W. Tichy and P. Marronetti, Phys. Rev. D 78, 081501 (2008).

[25] C. O. Lousto, H. Nakano, Y. Zlochower, and M. Campanelli, Phys. Rev. D 81, 084023 (2010); 82, 129902(E) (2010).

[26] A. E. Broderick, V. L. Fish, S. S. Doeleman, and A. Loeb, Astrophys. J. 697, 45 (2009).

[27] A. E. Broderick, V. L. Fish, S. S. Doeleman, and A. Loeb, Astrophys. J. 735, 110 (2011).

[28] V. I. Dokuchaev, Gen. Relativ. Gravit. 46, 1832 (2014).

[29] J. D. Schnittman, Astrophys. J. 806, 264 (2015).

[30] A. V. Maccio', G. Stinson, C. B. Brook, J. Wadsley, H. M. P. Couchman, S. Shen, B. K. Gibson, and T. Quinn, Astrophys. J. 744, L9 (2012).

[31] C. B. Brook, G. Stinson, B. K. Gibson, S. Shen, A. V. Macci, J. Wadsley, T. Quinn, and A. Obreja, Mon. Not. R. Astron. Soc. 443, 3809 (2014).
[32] M. Schaller, C. S. Frenk, R. G. Bower, T. Theuns, A. Jenkins, J. Schaye, R. A. Crain, M. Furlong, C. D. Vecchia, and I. G. McCarthy, Mon. Not. R. Astron. Soc. **451**, 1247 (2015).

[33] V. Springel et al., arXiv:1707.03397.

[34] P. Ullio, H. Zhao, and M. Kamionkowski, Phys. Rev. D **64**, 043504 (2001).

[35] D. Merritt, M. Milosavljevic, L. Verde, and R. Jimenez, Phys. Rev. Lett. **88**, 191301 (2002).

[36] D. Merritt, Phys. Rev. Lett. **92**, 201304 (2004).

[37] O. Y. Gnedin and J. R. Primack, Phys. Rev. Lett. **93**, 061302 (2004).

[38] G. Bertone and D. Merritt, Phys. Rev. D **72**, 103502 (2005).

[39] R. M. Buchholz, R. Schodel, and A. Eckart, Astron. Astrophys. **499**, 483 (2009).

[40] D. Merritt, *Dynamics and Evolution of Galactic Nuclei* (Princeton University Press, Princeton, NJ, 2013).

[41] R. Schdel, A. Feldmeier, N. Neumayer, L. Meyer, and S. Yelda, Classical Quantum Gravity **31**, 244007 (2014).

[42] E. Vasiliev, Phys. Rev. D **76**, 103532 (2007).

[43] S. L. Shapiro and J. Shelton, Phys. Rev. D **93**, 123510 (2016).

[44] C. Arina, T. Bringmann, J. Silk, and M. Vollmann, Phys. Rev. D **90**, 083506 (2014).

[45] V. Gammaldi, V. Avila-Reese, O. Valenzuela, and A. X. Gonzales-Morales, Phys. Rev. D **94**, 121301 (2016).

[46] C. S. Reynolds, Space Sci. Rev. **183**, 277 (2014).

[47] A. Belikov and J. Silk, Phys. Rev. D **89**, 043520 (2014).

[48] P. Young, Astrophys. J. **242**, 1232 (1980).

[49] G. D. Quinlan, L. Hernquist, and S. Sigurdsson, Astrophys. J. **440**, 554 (1995).

[50] R. P. van der Marel, Astron. J. **117**, 744 (1999).

[51] J. Binney and S. Tremaine, *Galactic Dynamics: Second Edition* (Princeton University Press, Princeton, NJ, 2008).

[52] A. G. Walker, Proc. Edinb. Math. Soc. **4**, 238 (1936).

[53] J. Bernstein, *Kinetic Theory in the Expanding Universe* (Cambridge University Press, Cambridge, U.K., 1988).

[54] E. D. Fackerell, Astrophys. J. **153**, 643 (1968).

[55] B. Carter, Phys. Rev. **174**, 1559 (1968).

[56] B. O'Neill, *The Geometry of Kerr Black Holes* (Courier Corporation, Chicago, IL, 2014).

[57] C. M. Will, *Classical Quantum Gravity* **29**, 217001 (2012).

[58] L. Sadeghian, Star Clusters and Dark Matter as Probes of the Spacetime Geometry of Massive Black Holes, Ph.D. thesis (2013), arXiv:1308.5378.

[59] T. L. Lee, T. Y. Li, and C. H. Tsai, Computing **83**, 109 (2008).

[60] G. P. Lepage, Report No. CLNS-80/447 (1980); G. P. Lepage, J. Comput. Phys. **27**, 192 (1978).

[61] W. H. Press and G. R. Farrar, Comput. Phys. **4**, 190 (1990).

[62] M. Galassi et al., *GNU Scientific Library Reference Manual* 3rd ed. (Network Theory Ltd., 2009). http://www.network-theory.co.uk/gsl/manual/

[63] L. Hernquist, Astrophys. J. **356**, 359 (1990).

[64] C. M. Will, Astrophys. J. **674**, L25 (2008).

[65] M. Fornasa and G. Bertone, Int. J. Mod. Phys. D **17**, 1125 (2008).

[66] E. Bulbul, M. Markevitch, A. Foster, R. K. Smith, M. Loewenstein, and S. W. Randall, Astrophys. J. **789**, 13 (2014).

[67] A. Boyarsky, O. Ruchayskiy, D. Iakubovskyi, and J. Franse, Phys. Rev. Lett. **113**, 251301 (2014).

[68] S. Dodelson and L. M. Widrow, Phys. Rev. Lett. **72**, 17 (1994).

[69] D. N. Spergel and P. J. Steinhardt, Phys. Rev. Lett. **84**, 3760 (2000).

[70] J. L. Feng, M. Kaplinghat, H. Tu, and H.-B. Yu, J. Cosmol. Astropart. Phys. **07** (2009) 004.

[71] M. Rocha, A. H. G. Peter, J. S. Bullock, M. Kaplinghat, S. Garrison-Kimmel, J. Onorbe, and L. A. Moustakas, Mon. Not. R. Astron. Soc. **430**, 81 (2013).

[72] S. L. Shapiro and V. Paschalidis, Phys. Rev. D **89**, 023506 (2014).

[73] A. X. Gonzalez-Morales, S. Profumo, and F. S. Queiroz, Phys. Rev. D **90**, 103508 (2014).

[74] F. Debbasch and W. A. van Leeuwen, Physica **388A**, 1818 (2009).

[75] I. S. Gradshteyn and I. M. Ryzhik, *Table of Integrals, Series and Products* (Academic Press, New York, 1980).