ANALYTIC STUDY OF MASS SEGREGATION AROUND A MASSIVE BLACK HOLE

A massive black hole (MBH) of mass M_• dominates the dynamics of stars within its radius of influence r_h = GM_•/σ² = 2.3(M_•/3 × 10⁶M_☉)(σ/75kms⁻¹)⁻²pc, where σ is the typical star velocity dispersion at r ≳ r_h and values quoted correspond to the Milky Way center (Alexander 2005). The steady-state distribution function (DF) of such stars, first derived (for simple, discrete stellar mass functions, MFs) by Bahcall & Wolf (1976, 1977; henceforth BW76,77), is useful in the study of galactic centers and possibly globular clusters. In the simple case where all stars have equal mass, BW76 showed that an n ∝ r^−7/4 density cusp forms; a somewhat flatter cusp was subsequently identified around SgA* (Alexander 1999; Genzel et al. 2003; Schödel et al. 2007).

Significant mass segregation is expected for realistic MFs (Miralda-Escudé & Gould 2000; Schödel et al. 2007; O'Leary et al. 2008, and references therein), as dynamical friction slows heavy stars that sink toward the MBH, while light stars are pushed outwards. This effect is essential in modeling various physical processes, such as tidal disruption (Lightman & Shapiro 1977; Magorrian & Tremaine 1999; Syer & Ulmer 1999) and gravitational wave emission (Hopman & Alexander 2006; Freitag et al. 2006; Hopman et al. 2007), and may be used for example to test if intermediate-mass black holes exist in globular clusters (Gill et al. 2008). Recently, Alexander & Hopman (2009, henceforth AH09) showed strong mass segregation in the limit of a strongly bimodal MF when the mass ratio is large and massive stars are rare.

In spite of extensive numerical studies of mass segregation around an MBH (AH09 and references therein), analytical modeling (BW77) is limited to simple, discrete, and extreme MFs where the most massive species alone constitutes more than half of the star population. Most results are limited to few stellar species and to a narrow mass range. In this Letter, we analytically study the steady-state single-star DF f within r < r_h, for an arbitrary, continuous MF. We confirm our results numerically for a wide range of MFs.

Following BW76,77, we assume (i) spherical spatial symmetry; (ii) isotropic velocities; (iii) Keplerian orbits; (iv) binaries are negligible; (v) small angle, uncorrelated, local scatterings; (vi) loss cone effects can be neglected; (vii) isothermal distribution of unbound stars with temperature μσ²; (viii) stars are destroyed when their specific energy drops below a threshold x ≡ −E/mσ² = x_max, with E the energy.

Consider the case where all stars have mass m. The dimensionless diffusion rate of stars through x is (BW76)

$$\begin{equation} Q(x) = \int _{-\infty }^{x_{\rm max}} \frac{f(x)f^{\prime }(x^{\prime })-f^{\prime }(x)f(x^{\prime })}{\left[\mbox{Max}(x,x^{\prime })\right]^{3/2}}\,dx^{\prime } = \mbox{const}\mbox{.} \end{equation} \tag{ 1 }$$

The boundary conditions are f(x < 0) = e^xm/μ and f(x>x_max) = 0. Equation (1) with these boundary conditions uniquely determines  f  and Q as shown below. The spatial number density of stars is given by (BW76)

$$\begin{equation} n(r) \propto \int _{-\infty }^{r_h/r} \,dx f(x)\sqrt{\frac{r_h}{r}-x} \mbox{.} \end{equation} \tag{ 2 }$$

Equation (1) describes the balance between diffusion toward the MBH and replenishment by new stars. BW76 found that Q is essentially determined by the diffusion rate at a bottleneck near x_max, where the replenishment rate diminishes. As Q ≃ 8/x_max is typically small, one can approximately solve Equation (1) by setting Q = 0.

We convert Equation (1) to an ordinary differential equation (ODE) by repeated operations of differentiation with respect to x and isolation of integral terms. This yields a nonlinear, fourth-order ODE for f,

$$\begin{equation} qf^2\left(B+4pA\right) = x^{3/2} Q C \mbox{.} \end{equation} \tag{ 3 }$$

Here we defined the local power-law index of f, p ≡ dln f/dln x = xw, the local power-law index of w ≡ (ln f)', q ≡ dln w/dln x, and the operators

$$\begin{eqnarray} A(f,x) & \equiv & 4q^2-5q-2 x^2w^{\prime \prime }/w\mbox{,}\\ B(f,x) & \equiv & 5q-12x^3w^{\prime \prime 2}/(w^{\prime }w)+8x^3w^{(3)}/w\nonumber \mbox{,}\\ C(f,x) & \equiv & \frac{pA^2 D}{3}\left[ \frac{B}{A^2} -\frac{4}{D} + \frac{3}{q}\right. \nonumber \\ & & \left. + \frac{2}{A}\left(10+2q-2p+\frac{5}{D} \right) \right] \nonumber \mbox{,} \end{eqnarray} \tag{ 4 }$$

where D(f, x) ≡ q + p − 5/2. The boundary conditions at x = 0 fix f and f' there, so the non-exponential solution f(x) is completely determined for given Q. A unique value of Q guarantees that f vanishes (for the first time) at x_max, proving the uniqueness of the steady-state solution (an "exercise left for the reader" by BW76).

For a power-law DF $f =f_0 x^{p_0}$, Equation (3) becomes

$$\begin{equation} 3f_0^2 \left(p_0-1/4\right) = Q x^{3/2-2p_0} p_0 \left(p_0-1 \right) \left(p_0-5/2 \right) \mbox{.} \end{equation} \tag{ 5 }$$

This compactly reproduces the BW76 results: although a finite, energy-independent flow requires p₀ = 3/4 (Peebles 1972, implying a flow of stars away from the MBH, see BW76), a steady-state with p₀ = 1/4 is closer to the actual DF because Q is negligibly small. Note that the power-law assumption becomes inconsistent at high energies, where $Q x^{3/2-2p_0}\propto x$ is large.

Henceforth, we assume Q = 0, justified for single mass stars if x_max ≫ 1. An exponential DF with q = 0 solves Equation (3), but does not satisfy the boundary conditions. For q ≠ 0, Equation (3) becomes a third-order nonlinear ODE,

$$\begin{equation} B+4pA=0 \mbox{,} \end{equation} \tag{ 6 }$$

with general solution of the form

$$\begin{equation} w=\frac{c_1}{x^{-3/2}-c_2} +c_3 + c_1 c_3 W(x;c_1,c_2,c_3) \mbox{,} \end{equation} \tag{ 7 }$$

with c₁, c₂, c₃ constants.

Useful results can be deduced directly from Equations (6) and (7), without determining the analytic properties of W. When p is much smaller (larger) than 1/4, Equation (6) becomes approximately B ≃ 0 (A ≃ 0), which corresponds to taking W (taking c₃) to zero in Equation (7). In this case f(x ≫ 1) ∝ exp [∫^xw(x')dx'] is either nearly constant, or (super-)exponential in x. In order to maintain reasonable stellar densities, there must be an approximate balance between the two terms, implying that p ≃ 1/4 far from the boundaries. Note that for p = 1/4, the six terms of A and B in Equation (4) precisely cancel in pairs.

For small x, where 0 ⩽ p = xw ≪ 1, Equation (6) becomes B ≃ 0; this corresponds to Equation (7) with W → 0, and c₃ = m/μ according to the boundary conditions at x = 0. At intermediate energies where approximately $f\propto x^{p_0}$, Equation (3) becomes −12w'^'2 + 8w⁽³⁾w' ≃ 0, with all other terms smaller by factors ≳5 (≳10) for |p₀| < 1/2 (|p₀| < 1/4), so w ≃ p₀(x + x₁)⁻¹ + x⁻¹₂. At high energies near the disruption energy, p is large and negative, Equation (6) becomes A ≃ 0 and so c₃ = 0; the boundary condition at x_max then yields c₂ = c₁/(3n) = x^−3/2_max with $n\in \mathbb {N}$. Smoothly combining these approximate solutions yields an approximation for the equal-mass case,

$$\begin{equation} f_{\rm eq} \simeq \left\lbrace\!\!\!\! \begin{array}{ll} \left(1+\frac{x m}{p_0 \mu } \right)^{p_0} e^{x/x_2} & \mbox{if $x<x_{ta}$;} \\ \begin{array}{ll} & \widetilde{f}\left[\frac{(\sqrt{x_{\rm max}}-\sqrt{x})^2}{x_{\rm max}+\sqrt{x_{\rm max} x} + x} \right]^n \exp \left[ \frac{3nx}{x_{\rm max}} \right. \\ & \quad\left. +\,\,2n\sqrt{3} \arctan \left(\frac{1}{\sqrt{3}} +\sqrt{\frac{4x}{3x_{\rm max}}} \right) \right] \end{array} & \mbox{if $x>x_{ta}$,} \end{array} \right. \end{equation} \tag{ 8 }$$

with constants x₂, x_ta and $\widetilde{f}$ determined, for example, by continuity of f, 3nx^−5/2_max/(x^−3/2 − x^−3/2_max) = p₀x⁻¹ + x⁻¹₂ and its first derivative evaluated at the turn-around energy x_ta. The results depend weakly on n; it is natural to fix n = 1. The previous paragraph suggests that p₀ → 1/4 as x_max increases. Numerically solving Equation (1) yields p₀ = (1/4)ρ(x_max), where the correction factor ρ ∼ 1 approaches unity as x_max becomes large, e.g., ρ(10⁴) = 1.18 and ρ(10⁸) = 1.08. It also depends weakly on the energy range used to measure p (here x^1/3_max < x < x^1/2_max, henceforth). Figure 1 illustrates the approximation.

Zoom In Zoom Out Reset image size

Consider stars with variable mass m in some range M_L < m < M_H. We generalize the discrete, multiple-mass version of Equation (1) (BW77) to the continuum limit,

$$\begin{eqnarray} Q(m) & = & \int _{M_L}^{M_H} m^{\prime }\,dm^{\prime } \int _{-\infty }^{x_{\rm max}} dx^{\prime } \, [\mbox{Max}(x,x^{\prime })]^{-3/2} \\ & & \times [m f(x,m)\partial _{x^{\prime }} f(x^{\prime },m^{\prime }) - m^{\prime } \partial _{x}f(x,m)f(x^{\prime },m^{\prime })] \mbox{,}\nonumber \end{eqnarray} \tag{ 9 }$$

where f is now the DF in x − m phase space. Assuming unbound stars with MF ξ(m), f(x < 0, m) = ξ(m)  e^xm/μ. The density n_m(r) is related to f (x, m) as in Equation (2), such that n_m(r ≳ r_h) ∼ ξ(m).

3.1. Negligible Flow

Consider the negligible flow limit Q → 0, which holds if the MF is not strongly dominated by light stars (AH09), as shown in Section 3.2. Dividing Equation (9) by mf(x, m) and differentiating with respect to m we find ∂_m[(mf)⁻¹∂_xf] = 0, implying that f must have the functional form

$$\begin{equation} f(x,m) = \xi (m)h(x)^m = \xi (m) e^{m k(x)} \mbox{,} \end{equation} \tag{ 10 }$$

with k ≡ ln h. Equation (10) implies that the local power-law index of f, p(x, m) ≡ dln f/dln x = mxk'(x), is linear in m for fixed x. This result, valid for negligible Q, generalizes the BW77 result p(x, m₁)/m₁ = p(x, m₂)/m₂ to a continuous MF.

Setting Q = 0 in Equation (9), using Equation (10), and repeatedly differentiating with respect to x and isolating integral terms, yields an ODE which, for q ≠ 0, becomes

$$\begin{equation} B+(p/P)A = 0 \mbox{.} \end{equation} \tag{ 11 }$$

Here, p, A, and B are functionals of f (m, x) and x, defined as in the equal-mass case (e.g., Equation (4)), and we defined

$$\begin{eqnarray} \hspace{-80pt} P(x,m) & \equiv & m/4M_0 \mbox{,} \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} M_L<M_0 & \equiv & \langle m^2 \rangle _f/\langle m \rangle _f = \langle m \rangle _{m f}<M_H \mbox{,} \end{eqnarray} \tag{ 13 }$$

with 〈X〉_Y ≡ ∫X(m)Y(m) dm/∫Y(m) dm averaging. The full distribution may be found by solving Equation (11) for k(x) under the boundary conditions k(x ⩽ 0) = x/μ and k → −∞ as x → x_max.

Note that P(x, M_L) ⩽ 1/4 and P(x, M_H) ⩾ 1/4. If the mass range is sufficiently narrow or ξ is sufficiently dominated by high masses such that P(x, M_H) ≃ 1/4, we recover the equal-mass case Equation (6) for the most massive stars. The full distribution then becomes

$$\begin{equation} f(x,m;\mu) \simeq \xi (m) f_{\rm eq}(x; \mu)^{m/M_H} \mbox{,} \end{equation} \tag{ 14 }$$

with f_eq(x; μ) the DF of equal-mass stars (Figure 1).

In the general case we may proceed as in Section 2: f will be constant, vanish, or diverge unless the two terms in Equation (11) are approximately balanced, i.e., p ≃ P. If P depends only weakly on x, this implies a power-law DF, f ∝ x^P. Thus, stars with mass equal to the weight-averaged mass M₀ tend to have a p = 1/4 cusp as in the equal-mass case. Taking some typical energy x^a_max with a ≲ 1 in Equation (13), we may calibrate a against numerical solutions of Equation (9). This yields a ≃ 1 for a wide range of MFs tested. Hence

$$\begin{equation} p(m) \simeq \rho (x_{\rm max})m/4M_0\mbox{,}\,\,\,\,\,\,\,\,\,M_0 = \langle m \rangle _{m \xi x_{\rm max}^{p(m)}} \mbox{,} \end{equation} \tag{ 15 }$$

and the full distribution is approximately

$$\begin{equation} f (x,m; \mu) \simeq \xi (m) f_{\rm eq}\left(x; \mu \right)^{m/M_0} \mbox{.} \end{equation} \tag{ 16 }$$

Equations (15) and (16) agree (to better than a factor of 2 in f) with numerical solutions of Equation (9) for various MFs, such as the discrete MFs tested by BW77 and AH09, the Salpeter MF (see Figure 2) and a wide range of power-law MFs (see Figure 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Equation (15) implies that for any MF, p(m) is a decreasing function of x_max. In particular, the power-law index of the most massive species, p_H, asymptotically approaches 1/4 as x_max → ∞, so f approaches Equation (14). In cases where we may approximate f ∼ ξ in Equation (13), for example if mξ is strongly peaked in a mass range where f/ξ = h^m varies little, we have M₀ ≃ 〈m〉_mξ.

For concreteness, consider unbound stars with a power-law MF ξ(m) ∝ m^α in some mass range M_L < m < M_H = ζM_L. In this case, we recover Equation (11) with

$$\begin{equation} P = -\frac{Mk}{4} \, \frac{\Gamma (2+\alpha,-M_Hk,-M_Lk)}{\Gamma (3+\alpha,-M_Hk,-M_Lk)} \mbox{,} \end{equation} \tag{ 17 }$$

and Γ(a, b, c) = ∫^c_bt^a−1e^−t dt the incomplete Γ-function. The power-law index p(m) may be found from Equation (15),

$$\begin{equation} \frac{\Gamma (3+\alpha,\frac{M_H}{m}\ln x_{\rm max}^{-ap}, \frac{M_L}{m}\ln x_{\rm max}^{-ap})}{\Gamma (2+\alpha,\frac{M_H}{m}\ln x_{\rm max}^{-ap}, \frac{M_L}{m}\ln x_{\rm max}^{-ap})} =\ln x_{\rm max}^{-a/4} \end{equation} \tag{ 18 }$$

(for one mass, then using p ∝ m). Indices p_L, p_H calculated from Equation (18) are illustrated in Figure 3 in the ζ − α plane. For small (very negative) α we may estimate p using the k → 0 limit where M₀ ∼ 〈m〉_mξ, whence⁶

$$\begin{equation} p_H \simeq \frac{\zeta (\alpha +3)(\zeta ^{2+\alpha }-1)}{4(\alpha +2)(\zeta ^{3+\alpha }-1)} \mbox{.} \end{equation} \tag{ 19 }$$

3.2. Non-Negligible Flow

When Q cannot be neglected, we may eliminate the x' integral in Equation (9) by repeated operations of differentiation with respect to x and isolation of integral terms. If we assume, in addition, a power-law energy dependence f ≃ ξ(m)x^p(m), we arrive at a generalization of Equation (5),

$$\begin{eqnarray} \xi (m) & {\int }_{M_L}^{M_H} & \left[ m\left(p^{\prime }+1\right)\left(p^{\prime }-\frac{1}{2}\right)+m^{\prime }p\left(\frac{3}{2}-p^{\prime }\right) \right] \nonumber \\ & & \, \times \, m^{\prime }\xi (m^{\prime }) x^{-\frac{3}{2}+p+p^{\prime }} dm^{\prime } \nonumber \\ & = & \frac{2}{3} p \left(p-1 \right) \left(p-\frac{5}{2}\right)Q(m) \mbox{,} \end{eqnarray} \tag{ 20 }$$

with abbrev. p = p(m) and p' = p(m'). For Q = 0 and p ∝ m, this reproduces the results of Section 3.1, p = m/4M₀.

Equation (20) suggests that p < 3/2 in all cases, otherwise the left-hand side becomes large and strongly x-dependent. Figure 3 shows that the p ≃ 3/2 limit is indeed realized if the MF is sufficiently broad and light-dominated. It corresponds to mass segregation even stronger than the saturation value p_H = 5/4 predicted by AH09 in the limit of a strongly bimodal, light-dominated MF where light stars are expected to assume p_L = 1/4.

A variant of Equation (20) with Q eliminated is

$$\begin{eqnarray} & {\int }_{M_L}^{M_H} & dm^{\prime } x^{p^{\prime }}m^{\prime }\xi (m^{\prime }) \left(\frac{3}{2}-p-p^{\prime }\right) \\ & & \times \left[ m\left(p^{\prime }+1\right)\left(p^{\prime }-\frac{1}{2}\right)+m^{\prime }p\left(\frac{3}{2}-p^{\prime }\right) \right] = 0 \nonumber \mbox{.} \end{eqnarray} \tag{ 21 }$$

Adopting a typical specific energy, x^b_max with b ≲ 1, one can solve this equation for p(m) with arbitrary ξ(m), an approximate procedure far simpler than solving Equation (9). However, Equation (21) typically becomes unstable when the underlying power-law assumption f ∼ x^p fails.

The linear scaling p(m) ∝ m remains approximately valid even when the flow is only marginally negligible. Using this Ansatz in Equation (21), p(m) becomes a root of the quadratic equation

$$\begin{equation} 4p^2 \langle m^2 \rangle - p(7-4p) m \langle m \rangle + \left(3/2-p\right)m^2=0 \mbox{,} \end{equation} \tag{ 22 }$$

with averages weighted by mξx^bp_max. For light-dominated MFs and massive stars with m ≫ {〈m〉, 〈m²〉^1/2}, the last term dominates and p(m) peaks at 3/2. For a strongly peaked MF where 〈m²〉 ≃ 〈m〉², Equation (22) yields two solutions: the Q = 0 limit p = m/4〈m〉, and a steep solution p = (3/2)(1 + 〈m〉/m)⁻¹.

Equation (22) suggests the following, numerically supported picture. For MFs peaked at relatively massive stars, the DF is approximately given by the Q = 0 result. For more light-dominated MFs, M₀ is shifted toward lower masses so the DF becomes gradually steeper at any given mass. If the MF is sufficiently extended (large ζ) and light-dominated (e.g., small α), each stellar species of mass m ≫ M_L eventually achieves maximal steepness p(m) ≃ 3/2. The transition between the Q = 0 and the saturation regimes involves in general an intermediate range of parameters in which the DF of this species is no longer a power law in x. We may roughly identify this region with the absence of real solutions to Equation (22),

$$\begin{equation} \delta (m)=(m/\langle m \rangle -5)^2-24(\langle m^2 \rangle /\langle m \rangle ^2-1) <0 \mbox{,} \end{equation} \tag{ 23 }$$

illustrated as areas delimited by solid curves in Figure 3. After reaching maximal steepness f(m) continues to evolve, but its maximal value remains ∼f(x = 1)x^3/2_max.

We have analyzed the distribution of an arbitrary, continuous MF ξ of stars around an MBH. For equal-mass stars, we derive a simple approximation for the DF f, see Equation (8) and Figure 1. For MFs not strongly dominated by light stars, where the mass flow can be neglected, we generalize the BW77 linear scaling p ∝ m to continuous MFs. We then derive approximate solutions for p (including its normalization, Equation (15)) and f (Equation (16)) and confirm them numerically; see Figures 2 and 3.

Our results provide a simple yet accurate alternative to solving the full integro-differential Equation (9). Equation (15) reproduces and generalizes previous Fokker–Planck calculations such as BW77, which were limited to a narrow range of parameters. As an illustration, consider SgA* with x_max = 10⁴ and a model MF with main-sequence stars, white dwarfs, neutron stars, and black holes, of masses M/M_☉ = 1, 0.6, 1.4, 10, and relative abundances 1:10⁻¹:10⁻²:10⁻³. Equation (15) then gives M₀ = 4.8M_☉ so p_H = 0.52 for black holes and p < 0.08 for the other species, in agreement with the numerical results of Hopman & Alexander (2006) and AH09.

The DF becomes gradually steeper with decreasing M₀, as long as δ in Equation (23) remains positive. For very light dominated MFs, the power-law assumption f ∼ ξx^p eventually fails at high masses, roughly where δ(m) < 0 (contour-enclosed regions in Figure 3). As the MF peak shifts to even lower masses, δ(m)>0 and the power-law behavior are eventually restored at the high-mass end and the DF remains very steep, peaking at p = 3/2. This behavior is illustrated in Figure 3.

Power-law DFs f ∝ x^p correspond to n(r) ∝ r^−3/2−p. As long as the MF is not very light dominated, stars of mass M₀ have p₀ = 1/4 (up to a ρ correction) and an equal-mass like cusp n₀ ∝ r^−7/4, while other stars have $n_m/n_0 \propto r^{-(m-M_0)/4M_0}$. If the MF is broad and heavy-dominated, light stars tend toward an energy-independent distribution n ∝ r^−3/2 while heavy stars approach n ∝ r^−7/4. In the opposite, maximally steep limit, heavy stars develop a cusp as steep as n ∝ r⁻³, harder than any reported previously, and the number of stars in the cusp depends (logarithmically) on the disruption cutoff.

We thank G. Van de Ven, D. Keres, and Y. Birnboim for helpful discussions. The research of U.K. was supported in part by NSF grant PHY-0503584, as a Friends of the Institute for Advanced Study member. The research of C.H. is supported by a Veni fellowship from the Netherlands Organization for Scientific Research (NWO). C.H. thanks the hospitality of the Institute for Advanced Study, where part of this research took place. T.A. is supported by ISF grant 928/06 and ERC grant 202996.