THE DISTRIBUTION AND ANNIHILATION OF DARK MATTER AROUND BLACK HOLES

Prompted by the recent paper by Bañados et al. (2009; BSW), there has been a great deal of interest in the potential of Kerr black holes to accelerate particles to ultra-relativistic energies and thus to probe a regime of physics otherwise inaccessible. The vast majority of this work has been analytic and thus largely limited to the most simple photon and particle trajectories in the equatorial plane. Here we present a more numerical approach that focuses on calculating the fully relativistic distribution function of massive test particles around a spinning black hole. With this distribution function and a simple model for the dark matter (DM) annihilation mechanism, we can then calculate the annihilation rate and observed spectrum as a function of black hole spin and observer inclination.

It has been noted repeatedly in recent works that the net energy gained through the Penrose process is quite modest, as is the fraction of collision products that might escape, and thus the astrophysical importance of the BSW effect is questionable (Jacobson & Sotiriou 2010; Bañados et al. 2011; Bejger et al. 2012; Harada et al. 2012; McWilliams 2013). We argue here that two primary factors (to our knowledge largely neglected in previous work) could greatly enhance the astrophysical relevance and observability of this annihilation. The first is an energy-dependent cross-section for DM annihilation. This could take many forms, the simplest of which is p-wave annihilation (Bertone et al. 2005; Chen & Zhou 2013; Ferrer & Hunter 2013), where the cross-section scales like the relative velocity, or a threshold energy, above which the cross-section increases greatly. This latter assumption is a natural choice for a model that includes multiple DM species, with the more massive particles intermediate products in the annihilation process toward gamma rays (see, e.g., Zurek 2014). Because gravity is the only known force capable of accelerating DM particles to high energies, it is possible that new annihilation channels could occur around black holes that are completely inaccessible everywhere else in the universe.

The other effect considered here is the relativistic enhancement of the density close to the black hole. This is due to the time dilation of observers near the horizon. In a steady-state system, one can think of dropping particles into the black hole from infinity at a constant rate ${{{\rm \Gamma }}}_{\infty }$ as measured by coordinate time t. To an observer near the black hole measuring proper time τ, an enhanced rate ${{{\rm \Gamma }}}_{\infty }({dt}/d\tau )$ is seen, with ${dt}/d\tau \gt 1$. For annihilation rates that scale like the density squared, the local annihilation rate will be enhanced by ${({dt}/d\tau )}^{2}$. Of course, the products will get redshifted on their way back out to an observer at infinity (McWilliams 2013), but we are still left with a net enhancement of ${dt}/d\tau$.

Even without this relativistic enhancement, numerous models also predict an astrophysical enhancement of the DM density in the galactic nucleus. Adiabatic growth of the central black hole will capture a large number of particles onto tightly bound orbits, growing a steep density spike as the black hole grows (Gondolo & Silk 1999; Sadeghian et al. 2013). Gravitational scattering off the dense nuclear star cluster will also lead to a DM spike (Gnedin & Primack 2004), similar to the classical two-body scattering result of Bahcall & Wolf (1976). At the same time, self-annihilation (Gondolo & Silk 1999) and elastic scattering (Fields et al. 2014; Shapiro & Paschalidis 2014) will act to flatten out this spike into a shallow core more similar to the unbound population.

Because our approach to this problem is predominantly numerical, we can easily include treatment of a range of black hole spins, particle distributions, and cross-sections, and not limit ourselves to special cases with analytic solutions. Therefore, we can calculate how often those extreme cases are likely to occur in a real astrophysical setting (for a notable exception to the analytic approaches of earlier work, see the exhaustive Monte Carlo calculations of Williams (1995, 2004) that explored the limits of the Penrose process in the context of Compton scattering and pair production in accretion disks and jets). Of particular interest has been the following question: for two particles, each of mass ${m}_{\chi }$, falling from rest at infinity and colliding near the black hole, what is the maximum achievable energy for an escaping photon? We find that this limit exceeds $12{m}_{\chi }$ for an extremal black hole with $a/M=1$, significantly higher than previously published values of $2.6{m}_{\chi }$ (Bejger et al. 2012). We explain the underlying reason for this discrepancy in a companion paper (Schnittman 2014).

2.1. Initial Conditions

The primary goal of this paper is to calculate the eight-dimensional phase-space distribution function ${df}({\bf{x}},{\bf{p}})$ of DM particles around a Kerr black hole. Two of these dimensions are immediately removed due to the assumption of a steady-state solution and stationarity of the metric, and the mass-shell constraint of the particle momentum, leaving us with ${df}(r,\theta ,\phi ,{p}_{r},{p}_{\theta },{p}_{\phi })$. This function is further reduced to five dimensions because axisymmetry removes the dependence on ϕ.

To calculate the distribution function, we first distinguish between two basic populations: the particles gravitationally bound and unbound to the black hole. The properties of the bound populations are more sensitive to underlying astrophysical assumptions, and will be discussed below in Section 2.4. The unbound population is more straightforward: we simply assume an isotropic, thermal distribution of velocities at a large distance from the black hole. Here, "large distance" is taken to be the influence radius r_infl of a supermassive black hole with mass M, and the DM velocity dispersion is set equal to the stellar velocity dispersion ${\sigma }_{0}$ of the bulge (thus the "unbound" population considered in this paper is still gravitationally bound to the galaxy, just not the black hole). From the "M–sigma" relation (Ferrarese & Merritt 2000) we take

$$\begin{eqnarray}&&M\approx 2\times {10}^{7}\;{M}_{\odot }{\left(\displaystyle \frac{{\sigma }_{0}}{100\;\mathrm{km}\;{{{\rm s}}}^{-1}}\right)}^{4}\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&{r}_{\mathrm{infl}}\equiv \displaystyle \frac{{GM}}{{\sigma }_{0}^{2}}\approx 8\;{{\rm }}\;{{\rm pc}}{\left(\displaystyle \frac{{\sigma }_{0}}{100\;\mathrm{km}\;{{{\rm s}}}^{-1}}\right)}^{2}.\end{eqnarray} \tag{ 2 }$$

In units of gravitational radii ${r}_{g}={GM}/{c}^{2}$, the influence radius is typically quite large: ${r}_{\mathrm{infl}}\approx {10}^{7}{M}_{7}^{-1/2}{r}_{g}$, where ${M}_{7}\;\equiv (M/{10}^{7}\;{M}_{\odot })$.

Given this outer boundary condition, we shoot test particles toward the black hole with initial velocities drawn from an isotropic thermal distribution with characteristic velocity ${\sigma }_{0}$. As we are only interested in the distribution function relatively close to the black hole, we can ignore any particle with an impact parameter greater than $\approx 1000\;{r}_{g}$. For those particles that we do follow, we calculate their geodesic trajectories with the Hamiltonian approach described in detail in Schnittman & Krolik (2013) and used in the radiation transport code Pandurata. A schematic of this procedure is shown in Figure 1. As the particle moves around the black hole and passes through different finite volume elements, the discretized distribution function ${df}({r}_{i},{\theta }_{j},{\bf{p}})$ is updated with appropriate weights.

Zoom In Zoom Out Reset image size

The great advantage of this Hamiltonian approach is that the integration variable is the coordinate time t in Boyer–Lindquist coordinates (Boyer & Lindquist 1967). Because this is the time measured by an observer at infinity, it determines the rate at which particles are injected into the system in the steady-state limit. Then the distribution function can be populated numerically by assigning a weight to each bin in phase space through which the test particle passes, with the weight proportional to the amount of time t spent in that volume. The process is repeated for many Monte Carlo test particles until the five-dimensional distribution function is completely populated.

2.2. Geodesics and Tetrads

Following Schnittman & Krolik (2013), we define local orthonormal observer frames, or tetrads, at each point in the computational volume. Depending on the population in question (i.e., bound versus unbound), it is convenient to use either the zero-angular-momentum observer (ZAMO; Bardeen et al. 1972) or the "free-falling from infinity observer" (FFIO) tetrads. In all cases we use Boyer–Lindquist coordinates (Boyer & Lindquist 1967), where the metric can be written

$$\begin{eqnarray}{g}_{\mu \nu }=\left(\begin{array}{cccc}-{\alpha }^{2}+{\omega }^{2}{\varpi }^{2} & 0 & 0 & -\omega {\varpi }^{2}\\ 0 & {\rho }^{2}/{{\rm \Delta }} & 0 & 0\\ 0 & 0 & {\rho }^{2} & 0\\ -\omega {\varpi }^{2} & 0 & 0 & {\varpi }^{2}\end{array}\right).\end{eqnarray} \tag{ 3 }$$

This allows for a relatively simple form for the inverse metric:

$$\begin{eqnarray}{g}^{\mu \nu }=\left(\begin{array}{cccc}-1/{\alpha }^{2} & 0 & 0 & -\omega /{\alpha }^{2}\\ 0 & {{\rm \Delta }}/{\rho }^{2} & 0 & 0\\ 0 & 0 & 1/{\rho }^{2} & 0\\ -\omega /{\alpha }^{2} & 0 & 0 & 1/{\varpi }^{2}-{\omega }^{2}/{\alpha }^{2}\end{array}\right),\end{eqnarray} \tag{ 4 }$$

with the following definitions:

$$\begin{eqnarray}&&{\rho }^{2}\equiv {r}^{2}+{a}^{2}{\mathrm{cos}}^{2}\theta ,\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{{\rm \Delta }}\equiv {r}^{2}-2{Mr}+{a}^{2},\end{eqnarray} \tag{ 5b }$$

$$\begin{eqnarray}&&{\alpha }^{2}\equiv \displaystyle \frac{{\rho }^{2}{{\rm \Delta }}}{{\rho }^{2}{{\rm \Delta }}+2{Mr}({a}^{2}+{r}^{2})},\end{eqnarray} \tag{ 5c }$$

$$\begin{eqnarray}&&\omega \equiv \displaystyle \frac{2{Mra}}{{\rho }^{2}{{\rm \Delta }}+2{Mr}({a}^{2}+{r}^{2})},\end{eqnarray} \tag{ 5d }$$

$$\begin{eqnarray}&&{\varpi }^{2}\equiv \left[\displaystyle \frac{{\rho }^{2}{{\rm \Delta }}+2{Mr}({a}^{2}+{r}^{2})}{{\rho }^{2}}\right]{\mathrm{sin}}^{2}\theta .\end{eqnarray} \tag{ 5e }$$

Unless explicitly included, we adopt units with $G=c=1$, so distances and times are often scaled by the black hole mass M.

The ZAMO tetrad can be constructed by

$$\begin{eqnarray}&&{{\bf{e}}}_{(\tilde{t})}=\displaystyle \frac{1}{\alpha }{{\bf{e}}}_{(t)}+\displaystyle \frac{\omega }{\alpha }{{\bf{e}}}_{(\phi )},\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{{\bf{e}}}_{(\tilde{r})}=\sqrt{\displaystyle \frac{{{\rm \Delta }}}{{\rho }^{2}}}{{\bf{e}}}_{(r)},\end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray}&&{{\bf{e}}}_{(\tilde{\theta })}=\sqrt{\displaystyle \frac{1}{{\rho }^{2}}}{{\bf{e}}}_{(\theta )},\end{eqnarray} \tag{ 6c }$$

$$\begin{eqnarray}&&{{\bf{e}}}_{(\tilde{\phi })}=\sqrt{\displaystyle \frac{1}{{\varphi }^{2}}}{{\bf{e}}}_{(\phi )},\end{eqnarray} \tag{ 6d }$$

where we designate tetrad basis vectors by $\tilde{\mu }$ indices, while coordinate bases have normal indices.

To construct the FFIO tetrad, the time-like basis vector ${{\bf{e}}}_{(\tilde{t})}$ is given by the 4-velocity ${u}^{\mu }={g}^{\mu \nu }{u}_{\nu }$ corresponding to a geodesic with ${u}_{t}=-1$, ${u}_{\theta }={u}_{\phi }=0$, and from normalization constraints,

$$\begin{eqnarray}&&{u}_{r}=-{\left[({\alpha }^{-2}-1)\displaystyle \frac{{\rho }^{2}}{{{\rm \Delta }}}\right]}^{1/2}.\end{eqnarray} \tag{ 7 }$$

Then the spatial basis vectors ${{\bf{e}}}_{(\tilde{i})}$ are constructed via a standard Gram-Schmidt method and aligned roughly parallel to the Boyer–Lindquist coordinate bases.

Any vector can be represented by its components in different tetrads via the relation

$$\begin{eqnarray}&&{\bf{u}}={{\bf{e}}}_{(\mu )}{{{\rm u}}}^{\mu }={{\bf{e}}}_{(\tilde{\mu })}{{{\rm u}}}^{\tilde{\mu }},\end{eqnarray} \tag{ 8 }$$

whereby the components are related by a linear transformation ${E}_{\tilde{\mu }}^{\mu }$:

$$\begin{eqnarray}&&{u}^{\mu }={E}_{\tilde{\mu }}^{\mu }{u}^{\tilde{\mu }},\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&{u}^{\tilde{\mu }}=[{E}^{-1}{]}_{\mu }^{\tilde{\mu }}{u}^{\mu }.\end{eqnarray} \tag{ 9b }$$

These ${u}^{\tilde{\mu }}$ are the components that we use for the tabulated distribution function.³ Because of the normalization constraints, we need only store three components of the 4-momentum in each spatial volume element, making the total dimensionality of the distribution function five: two space and three momentum.

In Pandurata, the geodesics are integrated with a variable time step fifth order Cash–Karp algorithm (Schnittman & Krolik 2013). This technique very naturally matches small time steps to regions of high curvature and thus areas of high resolution in the spatial grid. For each time step, a weight proportional to the coordinate time spent on that step is added to the distribution function for that particular volume of phase space. Because the particle typically remains within a single volume element for many time steps, we find that interpolation errors are small.

The spatial momentum components $\gamma {\beta }^{\tilde{i}}$ can be positive or negative and span many orders of magnitude. To adequately resolve the phase space and capture the relativistic effects immediately outside the black hole horizon, we find that on order ∼10³ bins are required in each dimension. If the entire phase-space volume were occupied, this would correspond to an unfeasible quantity of data. Fortunately, this volume is not evenly filled, so such a hypothetical five-dimensional array is in fact exceedingly sparse. In practice, we are able to use a dynamic memory allocation technique that only stores the non-zero elements of the distribution function. Yet even so, a well-resolved calculation can easily require multiple GB of data for a single distribution function, and to adequately sample this phase space requires on the order of ∼10⁹ test particles, with each geodesic sampled over thousands of time steps. Fortunately, this is a trivially parallelizable problem, so it is relatively simple to achieve sufficient resolution in a reasonable amount of time with a small computer cluster.

2.3. Unbound Particles

As mentioned above, for the unbound population, the outer boundary condition for the phase-space density at r_infl is relatively well-understood. The velocity distribution is thermal with characteristic speed ${\sigma }_{0}$.⁴ The spatial density of DM is measured from galactic rotation curves at kpc distances from the nucleus, and then must be extrapolated in to pc distances with a combination of observations and stellar profile modeling. For example, in the Milky Way the DM density near the Sun is 0.3 GeV cm⁻³, and the radial profile can be reasonably well-modeled with a simple $\rho \sim {R}^{-1}$ profile, giving a density of ∼10³ GeV cm⁻³ at r_infl. Inside of r_infl there is almost certainly an additional bound component to the DM distribution (Gondolo & Silk 1999), so the unbound population described here can best be understood as a strict lower bound on the phase-space density.

Outside of ∼100 ${r}_{g}$ the unbound population can be treated as a collisionless gas of accreting particles, as in Zeldovich & Novikov (1971). In the Newtonian limit, the density and velocity dispersion can be written as

$$\begin{eqnarray}&&n(r)={n}_{0}{\left(1+\displaystyle \frac{2{GM}}{{\sigma }_{0}^{2}r}\right)}^{1/2}\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&{\sigma }^{2}(r)={\sigma }_{0}^{2}\left(1+\displaystyle \frac{2{GM}}{{\sigma }_{0}^{2}r}\right).\end{eqnarray} \tag{ 11 }$$

In Figure 2 we show the spatial density of unbound particles as measured by a FFIO around a Kerr black hole with spin parameter $a/M=1$, as well as the mean particle momentum as measured in that frame. We find very close agreement to the Newtonian results all the way down to $r\sim 10{r}_{g}$. The deviation of the momentum from the Newtonian solution is due largely to the special relativistic terms proportional to the Lorentz boost γ.

Zoom In Zoom Out Reset image size

The proper density is governed by two competing relativistic effects. One is time dilation and the other is spatial curvature. Close to the black hole, the particle's proper time τ slows down relative to the coordinate time t measured by an observer at infinity, giving a large ${dt}/d\tau$. This has the effect of increasing the number density because, in a steady state, particles are injected into the system at a constant rate—as measured by an observer at infinity. The injection rate measured by an observer close to the black hole is higher by a factor of ${dt}/d\tau$, leading to her seeing a larger proper density.

In fact, the proper density would be even higher if it weren't for another important relativistic effect: the stretching of space around a black hole. Specifically, the Boyer–Lindquist radial coordinate element dr corresponds to a greater and greater proper distance as the observer approaches the horizon. This naturally gives a greater proper volume $d\tilde{V}$, shown as a solid curve in Figure 3. Again, we show the Newtonian value ${dV}/{dr}=4\pi {r}^{2}$ as a dashed curve. Because the particle interaction rates scale like ${n}^{2}\;v\;d\tilde{V}$, all these effects combine to increase the importance of reactions near the black hole.

Zoom In Zoom Out Reset image size

In Figure 4 we plot the momentum distributions of unbound DM particles, as observed by a FFIO in the equatorial plane, at a relatively large distance from the black hole: r = 100 M. Each one-dimensional distribution is calculated by integrating over the other two momentum dimensions. We also plot the momentum magnitude $\gamma | \beta |$ in panel (a). Because the particles all have relatively small velocities at infinity ${\beta }_{0}\approx {\sigma }_{0}/c\ll 1$, their velocities in the weakly relativistic region ${r}_{g}\ll r\ll {r}_{0}$ are given by $v\approx \sqrt{2{GM}/r}$, corresponding to $v\approx 0.14c$ for r = 100 M.

Zoom In Zoom Out Reset image size

For the three spatial components of the momentum distribution, we see a nearly isotropic velocity distribution with a few subtle but interesting deviations. First, we note how there is a slight deficit of particles with positive ${p}^{\tilde{r}}$. This is due to capture by the black hole of particles coming in from infinity with nearly radial trajectories. By definition, these particles also have small values of ${p}^{\tilde{\theta }}$ and ${p}^{\tilde{\phi }}$, depleting the distribution function in those dimensions around β = 0. While the distribution in the θ dimension is symmetric, note that the depletion in the ϕ distribution is offset to slightly negative values of ${p}^{\tilde{\phi }}$. This is due to the well-known preferential capture by Kerr black holes of retrograde particles with angular momenta aligned opposite to the black hole spin.

In Figure 5, we plot the phase-space distribution for the same boundary conditions as in Figure 4, but now at r = 2 M. The difference is quite dramatic, but all the features are essentially due to the same physical mechanisms. This close to the horizon, there is a very strong depletion of outgoing particles with ${p}^{\tilde{r}}\gt 0$, as most particles are captured by the black hole. The only particles that can avoid capture at this radius have prograde trajectories in the equatorial plane. Thus, the distribution is now peaked around ${p}^{\tilde{\theta }}=0$ instead of showing a deficit.

Zoom In Zoom Out Reset image size

There is also a strong peak near ${p}^{\tilde{\phi }}=1$ due to the relatively stable, long-lived prograde orbits that circle the black hole multiple times before getting captured or escaping back out to infinity. In fact, the distribution of coordinate momentum is significantly more lopsided to ${p}^{\phi }\gt 0$, but this is masked in Figure 5(d) because this distribution is measured by an observer with ${u}^{\phi }\gt 0$ herself. The sharp fall-off of the azimuthal distribution above ${p}^{\tilde{\phi }}\approx 1$ is due to the angular momentum barrier of the black hole. Particles with higher values of ${p}^{\tilde{\phi }}$ simply never reach this small radius.

To the best of our knowledge, these distribution functions have never been calculated before for a Kerr black hole. However, the particle number density can be determined analytically for a non-spinning Schwarzschild black hole, in the limit of ${\sigma }_{0}\ll c$. This allows at least one test of our numerical methods, although admittedly not a very strong one, as most of the interesting features are related to the far more complicated orbits around a spinning black hole. We follow the approach of Baushev (2009), who integrates the distribution function with fixed energy, carefully setting the angular momentum integration bounds based on which orbits are captured from a given radius. The results are shown in Figure 6, with our numerical calculation plotted as a red curve and the analytic result in black, showing perfect agreement. Note that Baushev's expression is given for a coordinate density rather than a proper density, which also explains the sharper peak at small r.

Zoom In Zoom Out Reset image size

2.4. Bound Particles

As mentioned above, the unbound population can be thought of as a lower limit on the total DM density. There will also likely be a substantial population of particles that are gravitationally bound to the black hole. As described in Gondolo & Silk (1999), the origin of the bound population is the adiabatic growth of the supermassive black hole on a timescale much longer than the typical orbital time. This physical mechanism can be understood as follows: as a marginally unbound DM particle passes within r_infl, a small amount of baryonic matter is accreted into this region, deepening the potential well just enough to capture the particle onto a marginally bound orbit. Once captured, the particle continues to orbit the black hole while conserving its orbital angular momentum as the black hole continues to gain mass. This has the effect of shrinking the radius of the orbit.

Over time, more particles are captured and subsequently migrate closer to the black hole, building up a steep density spike (Gondolo & Silk 1999). Inside of the inner-most stable circular orbit (ISCO), there is a sharp fall-off in the density spike due to plunging trajectories (Sadeghian et al. 2013). Here we do not attempt to solve for the slope of the density spike at large radii, but leave it as a free parameter, and fix the density at the influence radius as for the unbound population: ${n}_{\mathrm{bound}}(r)={n}_{0}{(r/{r}_{0})}^{-\alpha }$. Following Gondolo & Silk (1999), we also allow for the possibility of a density upper bound ${n}_{\mathrm{annih}}$ due to annihilation losses occurring over very long timescales.

To populate the phase-space distribution for the bound population, we follow a similar method as described above for the unbound particles, but instead of launching them from large radius with a limited range of impact parameters, now we launch them in situ with a isotropic thermal velocity distribution, as measured by a local ZAMO. These particles begin much closer to the black hole, so the relativistic Maxwell–Jüttner velocity distribution is used (Jüttner 1911), with the characteristic virial temperature ${{\rm \Theta }}(r)=1/2[1-{\epsilon }_{\mathrm{ZAMO}}(r)]$, where ${\epsilon }_{\mathrm{ZAMO}}(r)-1$ is the specific gravitational binding energy of the ZAMO.

Because many of the particles launched close to the black hole get captured, we first integrate their trajectories for a few orbital periods to ensure they are in fact on stable orbits. Only then do they contribute to the tabulated distribution function. Additionally, a small fraction of the test particles from the tail end of the velocity distribution will in fact be unbound, and these are similarly discarded.

As with the unbound distribution, for each step along its trajectory, the test particle contributes to the phase-space distribution a small weight proportional to the amount of time spent on that step. Yet now, instead of using the coordinate time dt, we use the proper time of the ZAMO frame from which the particles are launched, including an additional weight to ensure the appropriate radial form of the density distribution at larger radii.

In Figure 7 we plot the radial density distribution and mean relative momentum of the bound particles, as measured in the ZAMO frame, in the equatorial plane around a Kerr black hole with spin $a/M=1$. The density profile is constructed so that $\rho (r)\sim {r}^{-2}$ at large radii. We clearly see major differences relative to the unbound population shown in Figure 2. Because of the lack of stable orbits close to the black hole, the bound population declines inside $r\approx 4\;M$, which corresponds roughly to the mean radius of the ISCO for randomly inclined orbits around a maximally spinning black hole. This effect was described in detail for non-spinning black holes in Sadeghian et al. (2013). For equatorial circular orbits, only prograde trajectories are allowed inside of $r=9\;M$. This leads to all particles moving in roughly the same direction closer to the black hole, and explains why the relative momentum $\langle \gamma \beta {\rangle }_{\mathrm{rel}}$ does not increase nearly as fast for the bound population as it does for the unbound population, which allows plunging retrograde trajectories, and thus more "head-on" collisions.

Zoom In Zoom Out Reset image size

In Figure 8 we show the 2D density profile in the $x-z$ plane for both bound and unbound populations, for $a/M=0$ and $a/M=1$. The horizon in Boyer–Lindquist coordinates is plotted as a solid black line. For comparison purposes, the density scale is normalized to the mean value at r = 10 M. In reality, the density of the bound particles could be orders of magnitude greater at these radii (Gondolo & Silk 1999). The most obvious difference here is the depletion of bound orbits inside of the ISCO, which lies at $r=6\;M$ for non-spinning black holes. For spinning black holes, the radius of the ISCO is a function of the particle's inclination angle, ranging from $r=1\;M$ for prograde orbits in the equatorial plane, to $r=5.2\;M$ for polar orbits, and $r=9\;M$ for retrograde equatorial orbits.

Zoom In Zoom Out Reset image size

Inside of the ISCO, there is also the "marginally bound" radius, where particles with unity specific energy can exist on unstable circular orbits. This radius is also a function of inclination angle, and is plotted in Figure 8 as dotted curve. Inside of this orbit, no bound particles will be found (for improved visibility, we have left this region white, not black, as would be required by a strict adherence to the color scale). One interesting feature of Figure 8 is that the density of the unbound population around spinning black holes doesn't show any obvious θ-dependence. It appears that the enhanced density due to long-lived prograde orbits is almost exactly countered by the lack of retrograde orbits at the same latitude.

In Figure 9 we show the phase-space distribution for each of the momentum components, as measured by a ZAMO in the equatorial plane at large radius (r = 100M). Compared to the equivalent plot for the unbound distribution (Figure 4), we see a number of significant differences. First, the fact that these particles are bound requires that $E\lt 1$, and the imposed virial energy distribution results in mean velocities that are smaller than those of the unbound population by a factor of $\sim \sqrt{2}$. Second, because we require stable, long-lived orbits, there is a larger depletion around ${p}^{\tilde{\theta }}=0$ and ${p}^{\tilde{\phi }}=0$, as these trajectories are all captured by the black hole and thus do not contribute at all to the distribution function. Similarly, we see a larger asymmetry due to the preferential capture of retrograde orbits with ${p}^{\tilde{\phi }}\lt 0$.

Zoom In Zoom Out Reset image size

In Figure 10 we plot the same momentum distribution functions, now at r = 2 M. Here the contrast with the unbound population (Figure 5) is even greater. The only stable orbits at this radius are prograde, nearly circular, nearly equatorial orbits. This results in a relatively narrow distribution clustered around ${u}^{\tilde{\mu }}=[\sqrt{2},0,0,1]$ in the ZAMO frame. This narrower range in allowed velocities will have a profound impact on the shape of the annihilation spectrum, as we will see in the following section.

Zoom In Zoom Out Reset image size

Once we have populated the distribution function, we can calculate the annihilation rate given a simple particle-physics model for the DM cross-section. Again, it is simplest to work in the local tetrad frame. Including special relativistic corrections (Weaver 1976), the local reaction rate is given by the following:

$$\begin{eqnarray}&&R({\bf{x}})=\int {d}^{3}{{\bf{p}}}_{1}\int {d}^{3}{{\bf{p}}}_{2}\;f({\bf{x}},{{\bf{p}}}_{1})f({\bf{x}},{{\bf{p}}}_{2})\displaystyle \frac{{\gamma }_{\mathrm{rel}}}{{\gamma }_{1}{\gamma }_{2}}\;{\sigma }_{\chi }({\gamma }_{\mathrm{rel}}){v}_{\mathrm{rel}},\end{eqnarray} \tag{ 12 }$$

where ${\gamma }_{1}$ and ${\gamma }_{2}$ are the Lorentz factors of two particles as measured in the tetrad frame, ${v}_{\mathrm{rel}}$ is their relative velocity, and ${\sigma }_{\chi }$ is the annihilation cross-section (potentially a function of the relative velocity). $R({\bf{x}})$ has units of events per unit proper volume per unit proper time, so we multiply by $d\tau /{dt}$ to get the rate observed by a distant observer.

The distribution function $f({\bf{x}},{\bf{p}})$ is calculated numerically using the methods of Section 2. As discussed there, the numerical representation of f can have upwards of 10⁸ elements, so the direct integration of Equation (12) is generally not computationally feasible. Instead, we use a Monte Carlo sampling algorithm to pick random momenta for each particle with an appropriate weight based on the magnitude of f and the size of the discrete phase-space volume.

The spatial integration, however, is carried out directly, looping over coordinates r and θ. This is shown schematically in Figure 11. For each volume element, a large number (typically ∼ 10⁶) of pairs of particles are sampled, and for each pair, a center-of-mass tetrad is created. The total energy in the center-of-mass frame is given by

$$\begin{eqnarray}&&{E}_{\mathrm{com}}={m}_{\chi }\sqrt{2(1+{{\bf{u}}}_{1}\cdot {{\bf{u}}}_{2})},\end{eqnarray} \tag{ 13 }$$

where ${m}_{\chi }$ is the rest mass of the DM particle, and ${\bf{u}}={\bf{p}}/{m}_{\chi }$ is the particle 4-velocity.

Zoom In Zoom Out Reset image size

The 4-velocity of the center-of-mass frame is then given by

$$\begin{eqnarray}&&{{\bf{u}}}_{\mathrm{com}}=({{\bf{u}}}_{1}+{{\bf{u}}}_{2})/{E}_{\mathrm{com}}.\end{eqnarray} \tag{ 14 }$$

The center-of-mass tetrad is constructed with ${{\bf{e}}}_{(\tilde{t})}={{\bf{u}}}_{\mathrm{com}}$. The spatial basis vectors are totally arbitrary, as they are only needed to launch photons with an isotropic distribution in the center-of-mass frame. Two photons, labeled ${{\bf{k}}}_{3}$ and ${{\bf{k}}}_{4}$ in Figure 11, are launched in opposite directions, each with energy ${E}_{\mathrm{com}}/2$ in the center-of-mass frame. We then transform back to a coordinate basis for the geodesic integration of the photon trajectories to a distant observer.

As in Schnittman & Krolik (2013), for the photons that reach infinity, Pandurata can generate an image and spectrum of the emission region. An example in shown in Figure 12 for the annihilation signal from the unbound population around an extremal black hole, limiting the emission signal to the region $r\lt 100\;M$. While the flux clearly increases toward the center of the image, because the density and velocity profiles are relatively shallow (see Figure 2 above), the net flux is actually dominated by emission from large radii. These annihilation events are not very relativistic, so they produce a strong, narrow peak in the observed spectrum, centered at the DM rest mass energy.

Zoom In Zoom Out Reset image size

The annihilation events occurring closer to the horizon sample a much more energetic population of particles. Restricting ourselves to only those events where the center-of-mass energy is greater than 1.5× the combined rest mass of the annihilating particles, we can zoom in to the center of Figure 12. The result is shown in Figure 13, now focusing on the inner region within $r\lt 6\;M$. At these small radii, the effects of black hole spin become much more evident. One such effect is the characteristic shape of the Kerr shadow, defined by the impact parameter of critical photon orbits (Chandrasekhar 1983). The observed flux is clearly asymmetric, as the prograde photons originating from the left side of the image have a much greater chance of escaping the ergosphere and reaching a distant observer.

Zoom In Zoom Out Reset image size

There is another interesting feature of Figure 13 that we believe is novel to this work. Namely, the purple lobes emerging from the "mid-latitude" regions near the center of the image. These are regions of greater photon flux, albeit very highly redshifted. Recall, this image is created by considering only annihilations with a moderately high center-of-mass energy. Near the equatorial plane, extreme frame dragging ensures that the velocity dispersion is highly anisotropic, with most of the DM particles and their annihilation photons getting swept along on prograde, equatorial orbits. Above and below the plane, the DM distribution is more isotropic, leading to a more isotropic distribution of outgoing photons. Yet if one goes too far off the midplane, it becomes more difficult for the photons to escape. At the mid-latitudes, there is just enough frame dragging for photons to escape, yet not so much that they get deflected away from the observer.

The spectrum corresponding to this image is also plotted in Figure 13. Not surprisingly, the red and blue wings of the annihilation line shown in Figure 12 come from the most relativistic events. As pointed out by Piran & Shaham (1977), even reactions with very high center-of-mass energies will typically lead to photons with low energies as measured at infinity, thus explaining the red tail of the annihilation spectrum. The high-energy tail above $E=2{m}_{\chi }$ is due exclusively to Penrose-process reactions where one of the annihilation photons has negative energy and gets captured by the black hole (Penrose 1969; Piran et al. 1975).

Earlier analytic work predicted that the maximum energy attainable from the collisional Penrose process was $2.6{m}_{\chi }$ for particles falling from rest at infinity (Bejger et al. 2012; Harada et al. 2012). Because our calculation is fully numerical, it was able to reveal previously unknown trajectories leading to very high efficiencies with $E\gt 10{m}_{\chi }$, as seen in Figure 13. Closer inspection revealed that these high-energy photons are created when an infalling retrograde particle collides with a outgoing prograde particle that has just enough angular momentum to reflect off the centrifugal barrier, providing the necessary energy and momentum for the annihilation photon to escape the black hole (Schnittman 2014; Berti et al. 2015).

Due to the strong forward-beaming effects within the ergosphere, the escaping photon flux is highly anisotropic, with the peak flux and highest-energy photons emitted in the equatorial plane. Figure 14 shows the predicted annihilation spectra for observers at different inclination angles for the same DM profile as shown in Figure 13. Again, we restrict ourselves to the highest-energy reactions with ${E}_{\mathrm{com}}\gt 3{m}_{\chi }$.

Zoom In Zoom Out Reset image size

It is also instructive to plot the annihilation flux as a function of the emission radius. In Figure 15 we show both the observed flux (solid curves) and the flux that gets captured by the black hole (dashed curves) as a function of radius, integrated over all observing angles. The emission is further subdivided by the center-of-mass energy of the annihilating particles. Of course, the photons emitted closer to the black hole have a greater chance of getting captured. For the unbound population, the total escape fraction ranges from ${f}_{\mathrm{esc}}=93\%$ at r = 10 M down to ${f}_{\mathrm{esc}}(2\;M)=14\%$, and ${f}_{\mathrm{esc}}(1.1\;M)=0.25\%$. At small radii, these numbers are somewhat smaller than those calculated by Bañados et al. (2011), who only considered critical trajectories in the equatorial plane, where the escape probability is greatest. Yet at large radii, our distribution includes particles with typically greater impact parameters, and thus greater chance for escape.

Zoom In Zoom Out Reset image size

Another interesting feature of the curves in Figure 15 is the very sharp cut off above a critical radius for each energy bin. This is a natural consequence of conservation of energy. Because all unbound particles come in from rest at infinity with $E={m}_{\chi }$, the available kinetic energy in the center-of-mass frame is simply the gravitational potential energy ${{Mm}}_{\chi }/r$ at that radius. For example, to reach a center-of-mass energy of 10% above the rest mass energy, the particles must fall within $r\approx 10\;M$. Also note that inside $r\approx 4\;M$, most of the photons are captured, while outside of this radius, most escape. This is in close agreement with what we found for plunging orbits inside of the ISCO of a Schwarzschild accretion flow in Schnittman et al. (2013).

On the other hand, for the bound population of DM particles (Figure 15(b)), which by definition are not plunging, we find that the photon escape fraction is more than 90% at all radii, greatly increasing the relativistic effects observable from infinity. This is consistent with the classic calculation by Thorne (1974) which found that for thin accretion disks limited to circular, planar orbits outside the ISCO, the fraction of emission ultimately captured by the black hole was never more than a few percent, even for maximally spinning black holes where the majority of the flux emerges from extremely close to the horizon.

As we showed in Schnittman (2014), the peak energy attainable from particles falling in from infinity is a strong function of the black hole spin. Now, considering the full phase-space distribution function of the particles, we can see how the shape of the spectrum depends on spin. In Figure 16 we plot the flux seen by an equatorial observer, again limited to the high-energy annihilations with ${E}_{\mathrm{com}}\gt 3{m}_{\chi }$. For even marginally sub-extremal spins, the peak photon energy falls precipitously. As the spin decreases further, the number of collisions with ${E}_{\mathrm{com}}\gt 3{m}_{\chi }$ also decreases, thereby reducing the total flux observed. Lastly, the decreasing spin also increases the critical impact parameter for capturing prograde photons, making it harder for the annihilation flux to escape to infinity.

Zoom In Zoom Out Reset image size

Recall from Section 2.3 above that the density of the unbound distribution scales like $n\sim {r}^{-1/2}$. From the rate calculation in Equation (12) we see that the annihilation rate [events s⁻¹ cm⁻³] scales like $R(r)\sim {r}^{-3/2}$. Including the volume factor ${dV}=4\pi {r}^{2}{dr}$ we can write the differential annihilation rate as ${dR}/{dr}\sim {r}^{1/2}$. In other words, the unbound contribution to the annihilation signal diverges at large radius. In practice, the outer boundary can be set as the black hole's influence radius, typically 10^6–7r_g. This means that the observed signal will essentially be a delta function in energy, with only small perturbations from the relativistic contributions at small r, and thus measuring spin from annihilation lines would be a very challenging prospect indeed.

Two possible effects provide a way around this problem, each with its own additional uncertainties. One possibility is that the annihilation cross-section is a strong function of energy, increasing sharply above some threshold energy. This is admittedly rather speculative, and in conflict with leading DM models of self-annihilation (Bertone et al. 2005). On the other hand, we do not even know what the DM particle is, or if there are many DM species making up a rich "dark sector," with all the beauty and complexity of the standard model particles (Zurek 2014). One could easily imagine a DM analog of pion production via the collision of high-energy protons, in which case the only reactions could occur immediately surrounding a black hole, the ultimate gravitational particle accelerator. In this case, by construction the annihilation rate is dominated the region immediately surrounding the black hole.

Another possibility is that the DM density is dominated by a population of bound particles. As described above in Section 2.4, this population arises through the adiabatic growth of the black hole through accretion, capturing marginally unbound particles while also making the bound particles ever more tightly bound (Gondolo & Silk 1999; Sadeghian et al. 2013). This process will generally lead to a much steeper density profile, such as the n ∼ r⁻² distribution we use here. In this case, the differential reaction rate scales like ${dR}/{dr}\sim {r}^{-5/2}$ so the annihilation spectrum is now dominated by the particles at the smallest radii. In both cases—energy-dependent cross-sections and a large bound population—the relativistic effects described in Section 2.3 (expanded proper volume and time dilation) push the most important interaction region to even smaller radii, and thus the annihilation spectra are even more sensitive to the black hole spin.

In Figure 17 we show the annihilation spectra for both the bound and unbound populations for a variety of spins, now including emission out to $r=1000\;M$. The relative amplitudes are somewhat arbitrary, because we do not know what the relative densities of the two populations might be (see discussion below in Section 4), but it is almost certain that the bound population should dominate, possibly even by many orders of magnitude (Gondolo & Silk 1999). At the same time, the unbound signal will be even narrower and have a greater amplitude peak than shown here, as it is dominated by low-velocity particles at large radii. So while their overall amplitudes are uncertain, the detailed shapes of the spectra away from the central peak are relatively robust, depending only on the properties of geodesic orbits near the black hole.

Zoom In Zoom Out Reset image size

In this broad part of the spectrum, the bound and unbound signals show very different behavior. For non-spinning black holes, no particle can remain on a bound orbit inside of $r=4\;M$ (see Figure 8), so there are no annihilation photons coming from just outside the horizon, and these are the photons that produce the most strongly redshifted tail of the spectrum. As the spin increases and the ISCO moves to smaller and smaller radii, the line becomes steadily broader. On the other hand, the unbound particles are found all the way down to the horizon, where they can annihilate to highly redshifted photons regardless of the black hole spin.

Comparing Figures 5 and 10, we see that the unbound particles probe a much greater volume of momentum space at small radii. This in turn leads to a greater chance of producing the extreme Penrose particles that characterize the blue tail of the spectrum. Because all the bound particles are essentially on the same prograde, equatorial orbits, it is much more difficult to achieve annihilations with large center-of-mass energies, so the high-energy cut off in the spectrum is much closer to the classical result for a single particle decaying into two photons in the ergosphere (Wald 1974). In short, for bound particles the red tail of the spectrum is a better probe of black hole spin, while for the unbound population, the blue tail is the more sensitive feature. But in both cases, higher spin leads to a broader annihilation line.

In addition to the dependence on the DM density profile, the amplitude of the annihilation spectrum will also depend on the unknown DM mass and annihilation cross-section. At this point, it is only possible to use existing observations to set upper limits on these unknown parameters. One major obstacle that has plagued nearly all observational efforts to detect DM annihilation is the existence of more conventional astrophysical objects such as active galactic nuclei (AGNs), pulsars, and supernova remnants, all of which are powerful sources of high energy gamma rays. One solution to this problem is to focus on nearby dwarf galaxies, which are thought to have a high DM fraction and are not typically contaminated by AGN activity or significant star formation (Ackermann et al. 2011; note, however, the recent work by Gonzalez-Morales et al. (2014), which focuses on the contribution of black holes in dwarf galaxies).

Yet for our purposes, it turns out that the strongest upper limits actually come from the most massive galaxies with the most massive central black holes. Massive elliptical galaxies have the added advantage of being relatively quiescent both in nuclear activity and star formation (e.g., Schawinski et al. 2007). As mentioned above, the annihilation signal from the unbound population will be dominated by flux at large radii. It is difficult enough to spatially resolve even nearby black holes' influence radii with Hubble Space Telescope, much less gamma-ray telescopes, so any potential annihilation signal will tell us little about the black hole itself.

Prospects for detection of an unambiguous black hole signature improve if we consider annihilation models that include an energy dependence to the DM cross-section. For example, p-wave annihilation mechanisms will have cross-sections proportional to the relative velocity between the two annihilating particles (see Chen & Zhou 2013; Ferrer & Hunter 2013, and references therein). Unfortunately, from Equation (11) we see that this would only lead to an additional factor of ${r}^{-1/2}$ in the integrand of Equation (12), which would still be dominated by the contributions from large r.

This effect is shown in Figure 18, which plots the predicted spectra for two annihilation models: ${\sigma }_{\chi }(v)={{\rm const}}$ (black curves) and ${\sigma }_{\chi }(v)\propto v$ (red curves). The black hole spins considered are $a/M=0$ (dashed curves) and $a/M=1$ (solid curves), and in all cases only the unbound population is included. Integrating out to $r={10}^{4}\;M$, we see only a slight difference in the shape of the spectrum, with the ${\sigma }_{\chi }(v)\propto v$ model leading to a slightly broader peak (all curves are normalized to give a peak amplitude of unity).

Zoom In Zoom Out Reset image size

Another possible annihilation model is based on a resonant reaction at some energy above the DM rest mass, as suggested in Baushev (2009). If the cross-section increases sharply around a given center-of-mass energy, this would have the effect of focusing in on a relatively narrow volume of physical space around the black hole, as in Figure 15.

Alternatively, the cross-section could abruptly increase above a certain threshold energy, if new particles in the dark sector become energetically allowed, analogous to pion production via proton scattering. In either the resonant or threshold models for the annihilation cross-section, one might imagine a pair of heavier, intermediate dark particles getting created and then annihilating to two photons as in the direct annihilation model. If, for example, the mass of these intermediate particles is $1.5{m}_{\chi }$, then the observed spectrum would look like those plotted in Figures 13 and 16. With a significant increase in the cross-section above such an energy threshold, these relativistically-broadened spectra could in fact dominate over the narrow line component produced by the rest of the galaxy.

A less exotic option would be the simple density enhancement due to the bound population. If this is sufficiently large, it would easily dominate over the rest of the galaxy and also produce a characteristically broadened line sensitive to both black hole spin magnitude and orientation relative to the observer. Somewhat ironically, one of the things that could ultimately limit the strength of the annihilation signal from bound DM is annihilation itself. If the adiabatic black hole growth occurred at high redshift, then in the subsequent ∼10¹⁰ years, the bound population will get depleted via self-annihilation at an accelerated pace due to its high density (Gondolo & Silk 1999; Gonzalez-Morales et al. 2014).

On the other hand, if the black hole grows through mergers, or experiences even a single merger since the last extended accretion episode, it is quite likely that the bound DM population could get completely disrupted. The details of such an event are beyond the scope of this paper, but could be modeled by following test particles bound to each black hole through the merger, via post-Newtonian calculations (Schnittman 2010) or numerical relativity (van Meter et al. 2010).

The observational challenge is readily apparent: the black holes with the largest bound populations will tend to be in gas-rich galaxies with a lot of accretion and high-energy nuclear activity that could overwhelm the DM annihilation signal. The more massive black holes, residing in gas-poor quiescent galaxies, are also more likely to have lost their cloud of bound DM through a history of mergers. Even in the event that a gas-rich spiral galaxy hosts a quiescent nucleus, the black holes in those galaxies tend to have lower masses (Kormendy & Ho 2013).

While the relation between black hole mass and DM density is quite complicated for the bound population, it is relatively straightforward to calculate for the unbound population, which we can take as a lower bound on the DM density. Recall the influence radius r_infl is the distance within which the gravitational potential is dominated by the black hole, as opposed to the nuclear star cluster or DM halo. From Equation (2) we see that the influence volume scales like ${r}_{\mathrm{infl}}^{3}\sim {M}^{3/2}$, while the total mass enclosed is—by definition—on the order of M. If the DM and baryonic matter have similar profiles (by no means a certainty!), then more massive black holes should have lower surrounding DM density, with ${n}_{\mathrm{infl}\ }\sim {M}^{-1/2}$.

Because the unbound DM density falls off more rapidly outside the central core, the annihilation flux ${F}_{\mathrm{unbound}}$ will be dominated by the contribution from around r_infl, so we can estimate

$$\begin{eqnarray}&&{F}_{\mathrm{unbound}}\approx \displaystyle \frac{1}{{D}^{2}}{n}_{\mathrm{infl}}^{2}{\sigma }_{\chi }{v}_{\mathrm{infl}}{r}_{\mathrm{infl}}^{3}\sim {M}^{3/4}\;{D}^{-2},\end{eqnarray} \tag{ 15 }$$

with D the distance to the black hole and the mean velocity at the influence radius ${v}_{\mathrm{infl}}={\sigma }_{0}$.

If we consider a threshold energy annihilation model where all the flux comes from inside a critical radius ${r}_{\mathrm{crit}}\sim {{\rm few}}\times {r}_{g}$, then the density scales like ${n}_{\mathrm{crit}}\sim {n}_{\mathrm{infl}}{({r}_{\mathrm{infl}}/{r}_{\mathrm{crit}})}^{1/2}\sim {M}^{-3/4}$ while the relative velocity scales like ${v}_{\mathrm{crit}}\sim {\sigma }_{0}{({r}_{\mathrm{infl}}/{r}_{\mathrm{crit}})}^{1/2}\sim {M}^{0}$. The net flux then scales like

$$\begin{eqnarray}&&{F}_{\mathrm{unbound}}\approx \displaystyle \frac{1}{{D}^{2}}{n}_{\mathrm{crit}}^{2}{\sigma }_{\chi }{v}_{\mathrm{crit}}{r}_{\mathrm{crit}}^{3}\sim {M}^{3/2}\;{D}^{-2}.\end{eqnarray} \tag{ 16 }$$

In both cases, it appears that the brightest sources will be the closest, as opposed to the most massive.

Now consider the case where the annihilation signal is dominated by the bound contribution, the bound density is in turn limited by a self-annihilation ceiling as in Gondolo & Silk (1999), and there is a threshold energy above which the cross-section greatly increases. In this case, the flux is simply proportional to the total volume within the critical radius, so ${F}_{\mathrm{bound}}\sim {M}^{3}\;{D}^{-2}$. With this scaling, the greatest flux will actually come from more distant, more massive black holes. For example, NGC 1277, with a mass of $1.7\times {10}^{10}\;{M}_{\odot }$ and at a distance of 20 Mpc (van den Bosch et al. 2012), could give an observed flux over a thousand times greater than our own Sgr A*!

Recent works by Fields et al. (2014) and Gonzalez-Morales et al. (2014) have argued that current Fermi limits of gamma-ray flux from Sgr A* and nearby dwarf galaxies with massive black holes already place the strongest limits on annihilation from DM density spikes. Based on the arguments above, we believe that even stronger limits should come from more distant, massive galaxies. The other important advance presented in the present work is that, for either the energy-dependent cross-sections, or the steep density spikes, the annihilation signal will be dominated by the region closest to the black hole, and thus a fully numerical, relativistic rate calculation is absolutely essential.

Lastly, we should mention that gamma rays, while the primary observable feature explored in this work, are not the only promising annihilation product. High-energy neutrinos could also be produced in some annihilation channels, particularly those with energy-dependent cross-sections like p-wave annihilation (Bertone et al. 2005). While neutrinos obviously present many new detection challenges, the successful commissioning of new astronomical observatories like IceCube make this approach an exciting prospect (Aartsen et al. 2013). Furthermore, the non-DM backgrounds may contribute significantly less confusion in the neutrino sky.

As apparent in the previous section, there are still far too many unknown model parameters to allow for quantitative predictions of the annihilation flux from DM around black holes. Sadeghian et al. (2013) put it best: "there are uncertainties in all aspects of these models. However one thing is certain: if the central black hole Sgr A* is a rotating Kerr black hole and if general relativity is correct, its external geometry is precisely known. It therefore makes sense to make use of this certainty as much as possible." We have attempted to follow their advice to the best of our ability.

Thus, in order of decreasing confidence, the results in this paper can be summarized as follows.

• 1.  
  For a given DM density ${n}_{\mathrm{infl}}$ and velocity dispersion ${\sigma }_{0}$ at the black hole's influence radius, the fully relativistic, five-dimensional phase-space distribution has been calculated exactly for any black hole spin parameter, covering the region from r_infl all the way down to the horizon.
• 2.  
  Given this distribution function and a model for DM annihilation, the observed gamma-ray spectrum can be calculated by following photons from their creation until they are either captured by the black hole or reach the observer. Two important relativistic effects serve to increase the annihilation rate as compared to a purely Newtonian treatment: time dilation near the black hole effectively raises the density of the unbound population in a steady-state distribution being fed from infinity; and transforming from coordinate to proper distances greatly increases the interaction volume in the region immediately around the black hole (see Figure 3).
• 3.  
  Our numerical approach has unveiled previously overlooked orbits that can produce annihilation photons with extreme energies, far exceeding previous estimates for the maximum efficiency of the collisional Penrose process (Schnittman 2014). The peak energy attainable for escaping photons is a strong function of the black hole spin.
• 4.  
  The population of bound DM has also been calculated numerically, although this depends on two additional physical assumptions: a local isothermal velocity distribution with a virial-like temperature; and an overall radial power law for the density, as found in Gondolo & Silk (1999) and Sadeghian et al. (2013). Including only the long-lived stable orbits, we found that the density peaks in the equatorial plane somewhat outside of the ISCO, forming a thick, co-rotating torus around the black hole spin axis. Because the bound population is not plunging toward the horizon, the emerging flux has a much greater chance of escaping the black hole.
• 5.  
  The annihilation spectra from both the bound and unbound populations are sensitive to the spin parameter, but in opposite ways: the unbound spectrum varies mostly in the high-energy cut off, with higher spins allowing higher-energy annihilation products; the bound population moves closer and closer to the horizon with increasing spin, giving a stronger redshifted tail to the annihilation spectrum. Both bound and unbound spectra become more sensitive to observer inclination with increasing spin, as the spherical symmetry of the system is broken.
• 6.  
  For DM particle physics models with an energy-dependent cross-section (particularly one that increases with center-of-mass energy), the annihilation spectrum will be a more sensitive probe of the black hole properties. For DM models incorporating a rich population of dark sector species, black holes may be the most promising way to accelerate these particles and observe their interactions.
• 7.  
  The shape of the annihilation spectra is relatively robust, but the normalization is highly dependent on uncertain parameters such as the DM density profile and cross-section. If the unbound density profile follows the baryonic matter, with the shallow slopes seen in core galaxies, the observed flux should be a relatively weak function of black hole mass. If, on the other hand, the annihilation signal is produced by the most relativistic population within ${r}_{\mathrm{crit}}\sim {{\rm few}}\times {r}_{g}$, then the signal could scale like M³ and thus be dominated by the most massive black holes in the local universe.

While this paper has treated the bound and unbound particles separately, future work will also consider the self-interaction between these two populations (Fields et al. 2014; Shapiro & Paschalidis 2014), which may lead to a single, self-consistent steady-state distribution with density slope between −1/2 and −2. Future work will also focus on developing a robust framework in which we can use existing and future gamma ray observations to constrain various parameters of the particle physics (e.g., ${m}_{\chi }$, ${\sigma }_{\chi }(E)$, and the annihilation mechanism, i.e., line versus continuum) and astrophysical models (${n}_{\mathrm{infl}}$, the bound distribution normalization and slope, and the black hole mass, spin, and inclination). While initial work will focus on setting upper limits on reaction rates by looking at quiescent galaxies, our ultimate ambition is nothing short of an unambiguous detection of DM annihilation around supermassive black holes.

This work was partially supported by NASA grants ATP12-0139 and ATP13-0077. We thank Alessandra Buonanno, Francesc Ferrer, Ted Jacobson, Henric Krawczynski, Tzvi Piran, Laleh Sadeghian, and Joe Silk for helpful comments and discussion.