STELLAR ENERGY RELAXATION AROUND A MASSIVE BLACK HOLE

Observations of our own Galactic center (GC) reveal that stars there move on Keplerian orbits in the potential of a central compact object, which is consistent with a massive black hole (MBH) of mass M_• ∼ 4 × 10⁶ M_☉ (Eisenhauer et al. 2005; Ghez et al. 2005, 2008; Gillessen et al. 2009a, 2009b). It is widely assumed that the Milky Way is not unique in this sense, and that there is an MBH in the center of most, if not all, galaxies.

The GC, with its relatively low mass MBH, represents a sub-class of galactic nuclei where the dynamical relaxation time is expected to be shorter than the Hubble time (Alexander 2011), so that the stellar distribution near the MBH could have had time to reach a relaxed state that is independent of initial conditions. Moreover, because of its proximity, it is the only system to date where the dynamical state can be directly probed by observations. However, recent star counts in the GC (Buchholz et al. 2009; Do et al. 2009; Bartko et al. 2010) indicate that the radial density profile of the spectroscopically identified low-mass red giants on the ∼0.5 pc scale, thought to trace the long-lived population there, rises inward much less steeply (or even decreases) than is expected for a dynamically relaxed stellar population (Bahcall & Wolf 1976, 1977; Alexander & Hopman 2009; Preto & Amaro-Seoane 2010). The simplest interpretation is that, contrary to expectations, the GC is not dynamically relaxed by two-body interactions. This could be due, for example, to a cosmologically recent major merger event with another MBH, which carved out a central cavity in the stellar distribution (Merritt 2010), or due to the presence of faster competing dynamical processes that drain stars into the MBH (Madigan et al. 2011).

This interpretation is however further complicated by a recent study of the diffuse light density distribution in the GC (Yusef-Zadeh et al. 2012), which is thought to track the low-mass, long-lived main-sequence stars that are too faint and numerous to be individually resolved, and consist of the bulk of the stellar mass. In contrast with the stellar number counts, the diffuse light distribution appears to rise inward rather steeply down to the ∼O(0.01) pc scale, where destructive stellar collisions (not accounted for in the purely gravitational treatment of two-body relaxation) can plausibly suppress the stellar density profile (Alexander 1999).

These uncertainties about the dynamical state of the GC, and by implication, the state of galactic nuclei in general, have far-reaching ramifications. The possibility that the GC is unrelaxed severely limits the validity of extrapolating results derived from the unique observations of the GC to other galactic nuclei. While a relaxed system can be understood and modeled from first principles, independently of initial conditions, an unrelaxed one reflects its particular formation history, which is expected to vary substantially from galaxy to galaxy. This, and the possibility that the GC lacks the predicted high-density central stellar cusp, has in particular crucial implications for attempts to understand the dynamics of extragalactic gravitational wave (GW) sources and to predict their rates, since the low-mass MBH in the GC has long been considered the archetype of low-frequency extragalactic targets for planned space-borne GW observatories (Amaro-Seoane et al. 2007).

The study of dynamical relaxation around an MBH has a long history, beginning with the seminal studies of Peebles (1972), Bahcall & Wolf (1976, 1977), Cohn & Kulsrud (1978), and Shapiro & Marchant (1978), who used various simplifying approximations (Nelson & Tremaine 1999) to enable analytical and numerical approaches, and continuing more recently with studies of the time to approach steady state near an MBH in large N-body simulations (e.g., Preto et al. 2004; Preto & Amaro-Seoane 2010). A complementary approach of extracting the diffusion coefficients (DCs) from stellar energy evolution in N-body simulations was carried out in several studies (Lin & Tremaine 1980; Rauch & Tremaine 1996; Eilon et al. 2009; Madigan et al. 2011), leading to widely differing predictions about the dynamical state of galactic nuclei, especially very close to the MBH (Madigan et al. 2011). These open questions, and the fundamental importance of the process of orbital energy relaxation, motivate a re-examination of the analytic description of energy relaxation, coupled to detailed comparisons with carefully controlled simulations.

In this study we use analytical considerations and high-accuracy N-body simulations to address the following open questions: What is the relation between local energy relaxation and the approach to a global steady state? What are the properties of evolution due to energy relaxation on the short and intermediate timescales? How rare are extreme energy exchange events, and what is their contribution to relaxation? What is the physical meaning, and the correct value, of the Coulomb cutoffs introduced to control the formal divergences in the relaxation rates? In practical terms, these questions all reduce to the problem of relating the analytical/statistical description of local two-body interactions, to the global measurable properties of a stellar system. These issues must be addressed in order to better understand dynamical processes in galactic nuclei, improve their modeling, and reliably apply results from N-body simulations to real physical systems.

Our approach to these problems is complementary to that of Preto et al. (2004). Preto et al. (2004) simulated over a relaxation timescale a very large system with a non-Keplerian, and non-self-similar initial distribution function (DF) and potential, and a relatively low mass ratio Q = M_•/M_⋆. Preto et al. (2004) show that there is an agreement between the N-body and Fokker–Planck (FP) simulations by comparing the evolution of a small inner sub-system where the potential is dominated by that of the MBH and where the analytical/numerical solution of Bahcall & Wolf (1976; henceforth the BW cusp) is relevant.

Here, we carry out a suite of small N-body simulations of a system with realistically high values of Q and a potential that is dominated by the MBH, designed specifically to simplify the analytic treatment and to allow complete modeling of the time-dependent evolution on short timescales. This allows us to rigorously test and verify the theoretical predictions, and in particular the scaling of the relaxation with energy and with the slope of the cusp.

The choice of focusing on short timescales is to a large extent forced by the high computational cost of the simulations, since we model a high-mass-ratio system (Q = 10⁶, comparable to the mass ratio in the GC) with N ⩾ 10⁴ particles that are strongly bound to the MBH. The short-timescale behavior is however also of interest in itself, since dynamical processes near an MBH can be limited by the lifetime of the stars (for example, the massive stars in the GC), or by fast destructive processes (for example, collisional destruction or tidal disruption).

The elementary description of energy relaxation is by the energy transition probability, K_{E, J}(ΔE) (Equation (24)), which is the orbit-averaged probability per unit time, per unit energy, for the energy of a star to change from E to E + ΔE in a single scattering event. The resulting distribution of the accumulated change in energy over a finite time t, Δ_tE, is represented by the propagator, W_{E, J}(Δ_tE, t), which describes the energy distribution of an initially mono-energetic system of test stars, after a time t.

Past standard treatments of energy relaxation in stellar systems with the FP formalism (e.g., Bahcall & Wolf 1976) implicitly assumed that only the two lowest moments of the energy transition probability (i.e., the first and second DCs) play part in the evolution of the system. This is equivalent to assuming that the free propagator¹ evolves self-similarly as a Gaussian with width proportional to $\sqrt{t}$, that is, energy evolves by normal diffusion. However, as we show, the ∼1/|ΔE|³ divergence of the transition probability (Goodman 1983) implies that higher order DCs do play a role, which is dominant on short timescales, and diminishes with time, but still remains non-negligible even when the system relaxes to steady state. The implication is that energy relaxation progresses somewhat faster than $\sqrt{t}$, as a non-self-similar anomalous diffusion process.

Our approach of measuring energy relaxation in N-body simulations on short timescales forces us to develop a robust statistical model for describing these timescales. By treating the problem in Fourier space, we are able to derive the full propagator explicitly. It can then be evaluated either numerically in non-perturbative form, or analytically in a perturbative manner. This enables us to show that the propagator is a heavy-tailed DF that evolves as ${\sim} \sqrt{t\log (t)}$ on timescales much shorter than the FP energy diffusion timescale. This short-timescale behavior is an example of an anomalous diffusion process (Metzler & Klafter 2000), albeit a marginal one, since the anomality is only logarithmic. One consequence is the much higher probability of events that would be expected to be extremely rare for a Gaussian probability distribution. Another implication of the very slow convergence of this propagator to the central limit is that attempts to measure relaxation in relatively short N-body experiments can yield very misleading results, which reflect more the method used than the quantity probed. By analyzing the scaling properties of the energy propagator we suggest a robust statistical estimator, using fractional moments, which is independent of the high-energy cutoff. With it we measure the rate of energy relaxation on short timescales and confirm our analytical treatment, and in particular, the validity of the propagator and the DCs used in the standard FP approximation. Finally, we use the propagator to evolve a system by Monte Carlo (MC) simulations up to steady state, and thereby confirm the approximate validity of the FP results on long timescales.

The paper is organized as follows: In Section 2, we summarize the process of energy relaxation around a massive object in terms of a normal diffusion process described by the FP equation, and derive an analytical solution for the evolution of a stellar cusp to steady state. In Section 3, we present an elementary description of energy relaxation in terms of the transition probability and its integrated propagator, and describe the emergence of anomalous diffusion that ultimately approaches normal diffusion on long timescales. In Section 4, we formulate a method to analyze our N-body simulations and use it to confirm and calibrate our analytical description of relaxation. In Section 5, we briefly discuss the observable effects of anomalous diffusion on short-period stars orbiting the Galactic MBH, perturbed by a dark cusp of stellar remnants. We summarize and discuss our results in Section 6.

Determining the stellar distribution in the vicinity of an MBH is a classical problem in stellar dynamics, first discussed in the context of an MBH embedded in a globular cluster (Peebles 1972; Bahcall & Wolf 1976). The general time evolution of the phase-space DF of test particles with mass m_T, f_T(v, x, t), due to two-body interaction with field stars in a local velocity–mass distribution f_F(v, m_F), is given by the master equation (Goodman 1985)

$$\begin{eqnarray} \frac{D}{Dt}f_{T}(\mathbf {v},\mathbf {x},t) & = & \int d^{3}\mathbf {v}^{\prime }f_{T}(\mathbf {v}^{\prime },\mathbf {x},t)K(\mathbf {v}^{\prime },\mathbf {v})\nonumber \\ & & -f_{T}(\mathbf {v},\mathbf {x},t)\int d^{3}\mathbf {v}^{\prime }K(\mathbf {v},\mathbf {v}^{\prime }), \end{eqnarray} \tag{ 1 }$$

where D/Dt ≡ (∂/∂t) + v · (∂/∂x) − (∂ϕ/∂x) · (∂/∂v) and K(v, v') is the transition probability per unit time per unit velocity volume that a test star changed its velocity from v to v' = v + Δv as a result of a single encounter with a field star,

$$\begin{eqnarray} K(\mathbf {v},\mathbf {v}^{\prime }) & = & \frac{4G^{2}}{|\Delta v|^{5}}\int dm_{F}m_{F}^{2}\int _{\mathbf {w}\cdot \Delta \mathbf {v}}d^{2}\mathbf {w}\nonumber \\ & & \times f_{F}\left(\mathbf {w}+\mathbf {v}+\Delta \mathbf {v}\frac{m_{T}+m_{F}}{2m_{F}},m_{F}\right). \end{eqnarray} \tag{ 2 }$$

This description of stellar dynamics as a relaxation process assumes that changes in energy and angular momentum result from uncorrelated velocity deflections, which are local in time and space (e.g., Nelson & Tremaine 1999; Freitag 2008). Dynamical mechanisms that involve long-timescale correlations, e.g., resonant relaxation (Rauch & Tremaine 1996), are not taken into account here.

Assume that f_T and f_F are functions only of orbital energy E and time t; that the potential is dominated by a massive central object ϕ = GM_•/r; and that the field stars and the test stars have the same mass, M_⋆. The master equation (1) is then reduced to an equation of energy only (here and below the energy and potential of bound stars are defined as positive),

$$\begin{eqnarray} \frac{\partial N(E,t)}{\partial t} & = & \int d\Delta E\lbrace N(E-\Delta E)K_{E-\Delta E}(\Delta E)\nonumber \\ & & -N(E)K_{E}(\Delta E)\rbrace , \end{eqnarray} \tag{ 3 }$$

where K_E(Δ) is the orbit- and eccentricity-averaged energy transition probability, and N(E) is the stellar number density in energy space,

$$\begin{equation} N(E,t)=4\pi ^{2}J_{c}^{2}(E)P(E)f(E,t). \end{equation} \tag{ 4 }$$

By assuming weak encounters, ΔE ≪ E, and expanding N(E − ΔE) and K_{E − ΔE}(ΔE) around E − ΔE = E, the master equation can be expressed as a Kramers–Moyal expansion (e.g., Weinberg 2001)

$$\begin{equation} \frac{\partial N(E,t)}{\partial t}=\sum \frac{1}{n!}\left(-\frac{\partial }{\partial E}\right)^{n}N(E)\Gamma _{n}^{\Lambda }\left(E\right), \end{equation} \tag{ 5 }$$

where the DCs are the moments of the transition probability K_E(ΔE),

$$\begin{equation} \Gamma _{n}^{\Lambda }\left(E\right)=\int _{\Delta E_{\min }}^{\Delta E_{\max }}\Delta E^{n}K_{E}\left(\Delta E\right)d\Delta E, \end{equation} \tag{ 6 }$$

and where Λ = ΔE_max/ΔE_min. For n = 1, 2, Γ^Λ_n is proportional to the Coulomb logarithm log Λ (Equations (10) and (11)); for higher moments the dependence is much weaker, ∝(1 − Λ^{−⌊(n − 1)/2⌋}) (Equation (31)), and can usually be ignored.

By truncating the sum after the first two terms, the Kramers–Moyal expansion is reduced to the FP equation (e.g., Rosenbluth et al. 1957)

$$\begin{eqnarray} \frac{dN(E,t)}{dt} & = & \frac{1}{2}\frac{\partial ^{2}}{\partial E^{2}}\lbrace N(E,t)\Gamma _{2}^{\Lambda }(E)\rbrace \nonumber \\ & & -\frac{\partial }{\partial E}\lbrace N(E,t)\Gamma _{1}^{\Lambda }(E)\rbrace =-\frac{\partial F(E)}{\partial E}, \end{eqnarray} \tag{ 7 }$$

where F(E) is the stellar flux.

The FP approximation is valid only if the higher order terms can be neglected. This is usually justified by arguing that the n > 2 terms of the expansion are of the order of the n > 2 DCs themselves, which are smaller than the first two DCs by a factor of 1/log Λ (e.g., Binney & Tremaine 2008; Hénon 1960). However, as we show in Section 3.2, such a comparison is meaningful only when ∂/∂E ∼ 1/ΔE ∼ 1/E (Equation (31)), i.e., the change in energy has accumulated to an order unity change. This is a valid assumption only over long timescales; on short timescales the system cannot be described by the first two DCs only.

It is further commonly assumed that stars that are unbound to the MBH (E < 0) are drawn from a Maxwell–Boltzmann distribution,

$$\begin{equation} f(E)=\frac{n_{0}}{\left(2\pi \sigma _{0}^{2}\right)^{3/2}}e^{E/\sigma _{0}^{2}},\quad E<0, \end{equation} \tag{ 8 }$$

where n₀ and σ₀ are the stellar number density and velocity dispersion of this simplified model of the host galaxy. The inner energy boundary is determined by the maximal orbital energy E_d where stars can survive against any of the various disruption mechanisms that operate in the extreme environment near an MBH (e.g., tidal disruption, tidal heating, orbital decay by GW emission, stellar collisions),

$$\begin{equation} f(E)=0,\quad E>E_{d}. \end{equation} \tag{ 9 }$$

The steady-state solution (∂N/∂t = 0) of the FP equation can be simply obtained if a cusp solution f(E) = CE^p is assumed (e.g., Peebles 1972). Although this solution does not satisfy the boundary conditions (Equations (8) and (9)), it is in fact a reasonable approximation of the full solution far from the boundaries at 0 ≪ E ≪ E_d (e.g., Bahcall & Wolf 1976). For a full analytic solution see Keshet et al. (2009).

Approximating the DF by a pure power law, −1 < p < 1/2, the DCs are

$$\begin{equation} \Gamma _{1}^{\Lambda }(E)=-\frac{4(1-4p)\log \Lambda }{1+p}\frac{N(<E)}{Q^{2}P(E)}E\propto E^{p+1} \end{equation} \tag{ 10 }$$

and

$$\begin{equation} \Gamma _{2}^{\Lambda }(E)=\frac{16\left(3-2p\right)\log \Lambda }{(1+p)(1-2p)}\frac{N(<E)}{Q^{2}P(E)}E^{2}\propto E^{p+2}. \end{equation} \tag{ 11 }$$

The steady-state solution is then (Equation (7))

$$\begin{equation} (4p-1)(4p-3)\frac{2}{1-2p}=0, \end{equation} \tag{ 12 }$$

whose solutions are either p = 3/4 (Peebles 1972), which is unphysical since Γ^Λ₂(E) diverges with E_d → ∞ leading to a huge stellar flux away from the MBH (Bahcall & Wolf 1976), or p = 1/4, the so-called BW cusp solution (Bahcall & Wolf 1976). The BW cusp solution has zero net flow (F = 0) since the term in the equation that contains the second DC is independent of energy, and since there is no drift (Lightman & Shapiro 1977). This solution was confirmed by integrating the FP equation, both numerically (e.g., Bahcall & Wolf 1976; Cohn & Kulsrud 1978; Preto et al. 2004) and statistically (e.g., Shapiro & Marchant 1978), byMC simulations.

Stellar systems out of equilibrium evolve by redistributing their energy and angular momentum. When this proceeds by diffusion, the timescales associated with the various DCs generally reflect the dynamical evolution timescale. One such timescale is the energy diffusion time t_E,

$$\begin{equation} t_{E}(E)\equiv E^{2}/\Gamma _{2}^{\Lambda }(E), \end{equation} \tag{ 13 }$$

which expresses the timescale for a test star to change its orbital energy by order unity. Another commonly used quantity is the Chandrasekhar relaxation time (Chandrasekhar 1942; Spitzer 1987),

$$\begin{equation} t_{{\rm Ch}}\equiv \frac{\sigma ^{2}}{\langle \left.(\Delta v_{\parallel })^{2}\rangle \right|_{v=\sqrt{3}\sigma }}=0.34\frac{\sigma ^{3}}{G^{2}M_{\star }^{2}n\ln \Lambda }, \end{equation} \tag{ 14 }$$

where the second equality expresses the timescale for a typical test star with velocity $v=\sqrt{3}\sigma$, to change its kinetic energy by order unity, in a Maxwellian, isotropic, and homogeneous stellar system with stellar mass M_⋆, number density n, and velocity dispersion σ. The Coulomb logarithm is defined by the ratio of maximal and minimal impact parameters ln Λ = ln (p_max/p_min), where p_min = GM_⋆/σ² (small-angle approximation), and p_max is the size of the system.

These two timescales are related, but different, since t_Ch is local² in physical space (r, v), whereas t_E is local in (E, J) space, and hence is orbit-averaged over all (r, v) values along the orbit, and often samples a wide range of conditions in the stellar cluster. It is therefore to be expected that the values of the two timescales can be quite different, even though they express the same basic property of the system: Systems with short relaxation timescales evolve rapidly, and those with long timescales evolve slowly. For example, when the two timescale are compared at E ∼ σ², typically t_Ch ≫ t_E (Equation (21)).

The diffusive timescales cannot be used in themselves to estimate how long it takes an out-of-equilibrium system to approach steady state, assuming that such a state exists. This generally depends on the specific initial conditions, and the operative definition of convergence to steady state. Numeric experiments show that in many cases t_Ch, evaluated at the radius of MBH influence, is a reasonable estimator for the relaxation time to steady state, t_r.

We now proceed to use the FP equation to demonstrate analytically why that is so, and to derive the relation between t_E, t_Ch, and t_r in a power-law stellar cusp around an MBH. In Section 3.3, we verify that our analysis is not limited by the approximate nature of the FP treatment, by numerically integrating the system using MC simulations where the individual energy steps are generated using the isotropic propagator W_E(Δ_tE). Moreover, we show that strong encounters, which are ignored in the FP treatment, reduce the time to reach steady state by a factor of only ≲ 4/3.

Consider a background stellar population in steady state with a power-law cusp N_bg(E) = (5/4)N(E/E_p)^−9/4, where N is the total number of stars with E > E_h. The DCs (Equations (10) and (11)) are Γ^Λ₁(E) = 0 and Γ^Λ₂(E)∝E^9/4. Assume a small initial perturbation superimposed at energy E₀ on the steady-state background, N_p(E) = N_⋆δ(E − E₀). This system describes, for example, test stars that are scattered by a steady-state distribution of field stars. Over time, the perturbation will spread out until it relaxes to the ∝E^−9/4 steady-state distribution.

The dimensionless FP equation for the evolution of N_p is

$$\begin{eqnarray} \frac{dN_{p}(x,\tau)}{d\tau } & = & \frac{1}{2}\frac{\partial ^{2}}{\partial x^{2}}\lbrace N_{p}(x,\tau)x^{9/4}\rbrace , \end{eqnarray} \tag{ 15 }$$

where x = E/E_h and τ = t/t_E(E_h). We assume that the number of stars with E > E_h is constant (i.e., there is a reflective boundary at E_h). The perturbation then evolves as (see Equation (C11))

$$\begin{eqnarray} \eta (x,\tau) & \equiv \frac{N_{p}(x,\tau)}{N_{p}^{\infty }(x)}\!-\!1 & =\sum _{n=1}^{\infty }\! Q_{n}(x)Q_{n}(x_{0})e^{-j_{5,n}^{2}\tau /128}, \end{eqnarray} \tag{ 16 }$$

where N^∞_p(x) = N_bg(x)/N, $Q_{n}(x)=\sqrt{x/5}J_{4}(j_{5,n}x^{-1/8})/J_{4} (j_{5,n})$, J_n(x) is the Bessel function of the first kind, and j_{n, m} is the mth zero of J_n(x). Similar solutions for power-law background cusps with p > 0 are presented in Appendix C.

The time it takes for the small perturbation to decay, and for the system to converge back to its steady state, depends on the location of the initial perturbation, and on the operative definition of the convergence. The dynamical evolution back to steady state is slower for higher x₀ (closer to the MBH), in spite of the fact that the relaxation time decreases in that limit as t_E(x) ∼ x^−p (p = 1/4). This counterintuitive result can be understood by considering a perturbation of mass ΔM at x ⩾ x₀, where the enclosed steady-state mass scales as x^{p − 3/2}₀. It then follows the fractional perturbation scales as δM∝x^{3/2 − p}₀, and the rate for relaxing the perturbation scales as δM/t_E∝x^3/2₀. We therefore conservatively define the relaxation time as the time it takes a δ-function perturbation in the limit x₀ → ∞ to decay to an order unity perturbation, η(x_0,t_r) = 1, where

$$\begin{eqnarray} \eta (t) & = & \lim _{x_{0}\rightarrow \infty }\sum _{n=1}^{\infty }Q_{n}^{2}(x_{0})e^{-j_{5,n}^{2}t/(128t_{E})}\,\nonumber \\ & \approx & \frac{1}{1128\Gamma (5)^{2}}\frac{j_{5,1}^{8}}{J_{4}^{2}(j_{5,1})}e^{-j_{5,1}^{2}t/(128t_{E})}, \end{eqnarray} \tag{ 17 }$$

where here Γ is the Gamma function. This corresponds to the relaxation time

$$\begin{equation} t_{r}=\frac{128}{j_{5,1}^{2}}\log \left(\frac{1}{1128\Gamma (5)^{2}}\frac{j_{5,1}^{8}}{J_{4}^{2}(j_{5,1})}\right)t_{E}\simeq 11.1t_{E}. \end{equation} \tag{ 18 }$$

The conclusion that t_r ∼ O(10)t_E is robust, since it mostly reflects the time it takes the perturbation to propagate down to x = 1, rather than the exact details of the convergence criterion. For example, an alternative definition of t_r in terms of evaporation leads to a similar result. Consider the evaporation of a perturbation through the outer boundary x = 1, and define the relaxation time to steady state as the time it takes the number of stars still in the system to decrease by a factor P_⋆ = 0.1 (the conclusions depend only logarithmically on the exact value of P_⋆). In this setting the probability to find the test star in the system is given by (see Equation (C15))

$$\begin{equation} \frac{N_{p}(x>1,\tau)}{N_{\star }}=\frac{2}{\sqrt{x_{0}}}\sum _{n=1}^{\infty }\frac{J_{5}(j_{4,n})J_{4}(j_{4,n}/x_{0}^{1/8})}{j_{4,n}J_{5}^{2}(j_{4,n})}e^{-j_{4,n}^{2}\tau /128}, \end{equation} \tag{ 19 }$$

and the relaxation time (in the limit x₀ → ∞) is given by

$$\begin{equation} t_{r}\approx \frac{128}{j_{4,1}^{2}}\log \left(\frac{P_{\star }^{-1}j_{4,1}^{3}}{8J_{5}(j_{4,1})\Gamma (5)}\right)t_{E}\simeq 10t_{E}, \end{equation} \tag{ 20 }$$

similar to the relaxation time we obtained before.

This simple model can also explain the behavior seen in a variety of dynamical simulations of relaxation around an MBH (Bahcall & Wolf 1976, Figure 2; Hopman & Alexander 2006, Figure 1; Madigan et al. 2011, Figure 30): The initial relaxation to steady state is extremely rapid (super-exponential), while at later times it slows down to exponential decay. Equation (17) indeed shows that the perturbation converges to steady state as a series of exponential decaying terms with progressively faster rates 128/j²_{5, n} = {0.601, 1.19, 1.93, 2.81, 3.86, ...}. Thus, the higher terms vanish rapidly in a combined super-exponential rate, leaving last only the first and slowest term, which decays exponentially on a timescale ≳ t_E.

The relaxation time to steady state, t_r, can be related to the Chandrasekhar time t_Ch in the host galaxy (idealized here as an isotropic and homogeneous system). Beyond the MBH radius of influence at r_h = GM_•/σ²_∞, where σ_∞ is the asymptotic velocity dispersion, the stars are unbound to the MBH (E < 0), the velocity distribution is Maxwellian, the DF is given by Equation (8), and therefore t_Ch is given by Equation (14), with the asymptotic stellar density n₀. The typical orbital energy at r_h is ∼σ²_∞/2. Closer to the MBH, within the radius of influence where the dynamics are Keplerian, the steady-state DF for bound stars with E > σ²_∞ (∼r_h/2) is f ≈ 2f(0)(E/σ_∞)^p (e.g., Bahcall & Wolf 1976). The energy diffusion timescale is then (Equations (13) and (11))

$$\begin{equation} t_{E}\big(\sigma _{\infty }^{2}\big)=\frac{1}{64}\frac{Q^{2}}{N_{h}\log \Lambda }P_{h}\approx T_{{\rm Ch}}/44, \end{equation} \tag{ 21 }$$

where the orbital period P_h and star number N_h are evaluated at the energy scale of the MBH's sphere of influence, E_h = σ²_∞. It then follows that t_r ∼ 10t_E ∼ T_Ch/4 (Equations (18) and (21)). The large numeric prefactors that relate these different diffusion timescales highlight the point discussed above. While these timescales all generally reflect the tendency of the system to evolve, they measure different aspects of it: t_Ch is a local quantity in physical space, estimated outside the evolving region; t_E is an orbit-averaged quantity, estimated at some typical radius inside the evolving region; t_r describes the timescale to reach global steady state. In many cases of interest, t_r is not very different from the age of the stellar system (or the Hubble time, a commonly assumed upper limit). Inferences about the dynamical state of such systems should be mindful of the large differences between the meanings and values of the various diffusion timescales. The result t_r ∼ T_Ch/4 is consistent with numerical studies that show that a stellar system around an MBH approaches steady state in about T_Ch/3 to T_Ch/2 (e.g., Bahcall & Wolf 1976; Hopman & Alexander 2006).

The scaling of t_r with MBH mass M_• can be obtained from the relation t_r ≃ 10t_E, the empirical M/σ relation (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002),

$$\begin{equation} M_{\bullet }=M_{0}\left(\sigma /\sigma _{0}\right)^{\alpha }, \end{equation} \tag{ 22 }$$

and the stellar number parameterization N(> E) = Q(E/ξσ²)^−5/4, where ξ is an order unity factor. The time to steady state is then $t_{r}\!=\!(5\sqrt{2}\pi /64)\xi ^{-5/4}(M_{0}/M_{\bullet })^{3/\alpha }(GM_{\bullet }/\sigma _{0}^{3}) Q/\log Q$. For α = 4, M₀ = 1.3 × 10⁸ M_☉, and σ₀ = 200 km s⁻¹ (Tremaine et al. 2002),

$$\begin{equation} t_{r}=(7.8\times 10^{9}\,\mathrm{yr})\xi ^{-5/4}\left(\frac{M_{\star }}{M_{\odot }}\right)^{1/4}\frac{(Q/10^{7}){}^{5/4}}{(\log Q/\log 10^{7})}. \end{equation} \tag{ 23 }$$

This indicates that stellar cusps around lower-mass MBHs (M_• ≲ 10⁷ M_☉ for M_⋆ = 1 M_☉ and ξ = 1), which have evolved passively over a Hubble time, should be dynamically relaxed.

This estimate of the relaxation time for small perturbations is broadly consistent with the factor 1–2 longer one found for major perturbations by Preto et al. (2004), as well as the results of Antonini et al. (2012), if the artificial slowing of relaxation due to numeric force smoothing in N-body simulations is taken into account.

The transition probability encodes the full description of energy evolution in the system due to two-body encounters on all energy scales. For hyperbolic Keplerian interactions between two stars, its general form is K_{r, E} ∼ F_{r, E}(ΔE)/|ΔE|³, where F_{r, E} is a non-diverging function of ΔE, the energy change in an encounter (see Appendix A). The singularity of K_{r, E} at ΔE → 0 is never reached in practice, because the finite size of the system sets a lower limit on ΔE. In fact, an even stronger limit is implied by the restriction to hyperbolic orbits. A test star on an orbit with energy E has radial orbital period P(E) ∼ GM_•E^−3/2 and a typical velocity $v\sim \sqrt{E}$. An encounter with another star at impact parameter b leads to an exchange of energy of order ΔE ∼ GM_⋆/b and lasts a time ∼b/v. Since the weakest encounter that can still be approximated by a hyperbolic one, must have b/v < P, it follows that ΔE_min ∼ E/Q (i.e., Δ_min ∼ 1/Q for Δ ≡ ΔE/E) (Bahcall & Wolf 1976). This approximate effective low-ΔE cutoff assumes that the contribution of distant and slow interactions is negligible. It is estimated below empirically by N-body simulations (see Section 4.2).

Kinematics introduce an inherent high-ΔE limit (Goodman 1985). Strong encounters occur in the impulsive limit, where the distance between the test star at r and the background star at r_b is very small. The energy is then exchanged over a short time and the test star's position remains almost fixed. In that limit, the energy one star can transfer to the other is at most its total kinetic energy. Let the test star be on an initial Keplerian orbit with semimajor axis a, energy E = GM_•/2a, eccentricity e, and let its eccentric anomaly at the moment of the encounter be ω, so that the radius then is r = a(1 − ecos ω). The most bound orbit it can be deflected to at fixed r has a'_min = r/2 and E'_max > E. Thus, the maximal change in the relative energy of the test star is Δ^{(+ )}_max = max |E' − E|/E = (2a/r − 1). Conversely, the least bound orbit the test star can be deflected to at fixed r involves kicking a background star with an initial energy E_b to a maximally bound orbit with a'_b = r_b/2 ≈ r/2. Since typically the least bound (or unbound) background stars have E_b ∼ 0, the maximal change in the relative energy in this case is Δ^{(− )}_max ≈ 2a/r. Figure 1 shows the full transition probability for a specific value of a/r and a power-law DF; the kinematic high-ΔE limit is clearly seen. The high-ΔE limit is ultimately set by physical (inelastic) collisions between stars of finite size, where orbital energy is not conserved due to tidal interactions, mass loss, or stellar destruction. Orbit-averaging over the phase of the test star yields 〈Δ^{(+ )}_max〉 = 1 and 〈Δ^{(− )}_max〉 = 2, close to the commonly assumed value Δ_max = 1.

Zoom In Zoom Out Reset image size

The kinematic limits have several noteworthy properties. (1) The averaged limits are consistent with the small-angle regime, conventionally defined by an impact parameter b > 2GM_⋆/〈v〉² (e.g., Binney & Tremaine 2008), which corresponds to Δ < 1. (2) In the limit e → 1, Δ diverges, and therefore the eccentricity-averaged transition probability has finite contributions from arbitrarily large Δ. In particular, the isotropic average transition probability, and the resulting isotropic propagator, involves integrations over Δ → ∞. (3) The actual value of Δ_max may be lower than the kinematical one, when destructive physical collisions between stars play a role. In that case the maximal energy exchange is given by Δ^⋆_max ∼ 4a/QR_⋆ ≈ 45(Q/4 × 10⁶)⁻¹(a/pc)(R_⋆/R_☉)⁻¹, where R_⋆ is the stellar radius. However, this is relevant only extremely close to the MBH. For example, in the GC Δ_max < 1 for a ≲ 22 mpc, and even then the effect of the collisional limit on the relaxation of the surviving stars is small; for example at a ∼ 0.22 mpc, log Λ is reduced by only a factor of two.

In order to measure short-timescale evolution in N-body simulation, it is necessary to derive the energy propagator in the short time limit. Consider a test star with an initial energy E and angular momentum J, where the energy and angular momentum change as a result of two-body interactions over a period of time P(E) < t < t_E. In the small deflection angle limit (equivalently, ΔE/E < 1), the orbit-averaged energy transition probability, K_{E, J}(ΔE) = ∫dΔJK_{E, J}(ΔE, ΔJ), can be approximated as (see Appendix A and Goodman 1983)

$$\begin{equation} K_{E,J}(\Delta E)=\frac{\Gamma _{2}(E,J)+\Gamma _{1}(E,J)\Delta E+O(\Delta E^{2})}{2\left|\Delta E\right|^{3}}, \end{equation} \tag{ 24 }$$

for ΔE_min < |ΔE| < E, where Γ₂(E, J) and Γ₁(E, J) are given by Equations (A6) and (A7). Note that Γ_{1, 2} are shape parameters of the transition probability, K_{E, J}(ΔE), and are not by themselves the DCs Γ^Λ_{1, 2}, which are related to them by

$$\begin{equation} \Gamma _{1,2}^{\Lambda }\left(E,J\right)=\Gamma _{1,2}(E,J)\log \Lambda. \end{equation} \tag{ 25 }$$

This distinction is important because log Λ ∼ log Q is typically a large, O(10) factor, which is determined by our assumptions that strong encounters are ignored. We address this assumption and estimate the contribution of strong encounters in Section 3.3 below. However, we show here that the short-timescale evolution of the system does not depend on the high-ΔE limit.

Assume that t is short enough so that the background remains fixed, and changes in E and J are small enough so that K_{E, J}(ΔE) is approximately constant along the trajectory of the test star in (E, J) space. Let W_{E, J}(Δ_tE, t) be the energy propagator. In that case the master equation is

$$\begin{eqnarray} \frac{\partial W_{E,J}(\Delta _{t}E,t)}{\partial t} & = & \int K_{E,J}(\Delta E)W_{E,J}(\Delta _{t}E-\Delta E,t)d\Delta E\nonumber \\ & & -q_{2b}(E,J)W_{E,J}(\Delta _{t}E,t), \end{eqnarray} \tag{ 26 }$$

where W_{E, J}(Δ_tE, 0) = δ(Δ_tE) and q_2b is the total rate at which two-body encounters scatter the test star,

$$\begin{equation} q_{2b}\left(E,J\right)=\int d\Delta EK_{E,J}(\Delta E)=\frac{1-\Lambda ^{-2}}{2\Delta E_{\min }^{2}}\Gamma _{2}(E,J). \end{equation} \tag{ 27 }$$

The energy propagator W_{E, J}(Δ_tE, t) can then be used to integrate the master equation in an MC scheme (e.g., Shapiro & Marchant 1978). Note that the Λ⁻² ∼ Q⁻² term in Equation (27) expresses the very low probability of the neglected large-angle scattering events when the transition probability is approximated as a power law (Equation (24)). It can serve as an order-of-magnitude estimate of the contribution of such of rare events to the total encounter rate in a real system.

The Fourier transform of Equation (26) is

$$\begin{equation} \left(\frac{\partial }{\partial t}-\tilde{K}_{E,J}(k)+q_{2b}(E,J)\right)\tilde{W}_{E,J}(k,t)=0, \end{equation} \tag{ 28 }$$

where the solution is simply

$$\begin{eqnarray} W_{E,J}(\Delta _{t}E,t) & = & \int _{-\infty }^{\infty }\frac{dk}{2\pi }e^{ik\Delta _{t}E}e^{(\tilde{K}_{E,J}(k)-q_{2b})t} \end{eqnarray} \tag{ 29 }$$

and $\tilde{K}_{E,J}(k)$ is the Fourier transform of K_{E, J}(ΔE),

$$\begin{eqnarray} \tilde{K}_{J,E}(k) & = & q_{2b}-ik\Gamma _{1}^{\Lambda }-\frac{k^{2}}{2}\Gamma _{2}^{\Lambda }+\sum _{n=3}^{\infty }\frac{(-ik)^{n}}{n!}\Gamma _{n}^{\Lambda }, \end{eqnarray} \tag{ 30 }$$

where Γ^Λ_n are the nth (non-central) moments of K_{J, E}(ΔE). The higher odd and even moments are

$$\begin{equation} \Gamma _{2n+1}^{\Lambda }=\frac{E^{2n}\Gamma _{1}}{2n}\frac{1-\Lambda ^{-2n}}{\Delta \epsilon _{\max }^{2n}},\quad \Gamma _{2n+2}^{\Lambda }=\frac{E^{2n}\Gamma _{2}}{2n}\frac{1-\Lambda ^{-2n}}{\Delta \epsilon _{\max }^{2n}}. \end{equation} \tag{ 31 }$$

Note that Equation (30) has an analytical expression in terms of cosine and sine integrals (see Equation (B2)). Thus, a full non-perturbative evaluation of the propagator W_{E, J}(ΔE, t) is obtained by numerically evaluating Equation (29). This is equivalent to including in the master equation DCs of all orders (i.e., all moments of the transition probability), although calculated using an approximate transition probability, which incorporates the large-angle limit only as an effective cutoff. This is to be contrasted with the standard FP treatment where only the first two DCs are included (often also only as approximations). For example, this truncation is valid for a Gaussian transition probability, where the higher order terms fall off as Λ⁻²ⁿ instead of (1 − Λ⁻²ⁿ)/Δ²ⁿ_max.

Neglecting the contribution of the n > 2 moments in Equation (30) results in the FP propagator

$$\begin{equation} W_{E,J}^{{\rm FP}}(\Delta _{t}E,t)=\frac{1}{\sqrt{2\pi \Gamma _{2}^{\Lambda }t}}e^{-\frac{(\Delta _{t}E-\Gamma _{1}^{\Lambda }t)^{2}}{2\Gamma _{2}^{\Lambda }t}}, \end{equation} \tag{ 32 }$$

which can be related to the non-perturbative energy propagator (Equation (29)) via Equation (30),

$$\begin{eqnarray} &&W_{E,J}(\Delta _{t}E,t) = W_{E,J}^{{\rm FP}}(\Delta _{t}E,t)\nonumber \\ &&\quad\quad\times \Bigg\lbrace 1+-\frac{1}{2}\frac{1}{3!}\frac{1}{\log \Lambda }\frac{E\Gamma _{1}}{\Gamma _{2}}\sqrt{\frac{t_{E}}{t}}H\! e_{3}\left(\frac{\Delta _{t}E-t\Gamma _{1}^{\Lambda }}{\sqrt{t\Gamma _{2}^{\Lambda }}}\right)\nonumber \\ &&\quad\quad +\frac{1}{2}\frac{1}{4!}\frac{1}{\log \Lambda }\frac{t_{E}}{t}H\! e_{4}\left(\frac{\Delta _{t}E-t\Gamma _{1}^{\Lambda }}{\sqrt{t\Gamma _{2}^{\Lambda }}}\right)+\dots \Bigg\rbrace, \nonumber \\ & & \end{eqnarray} \tag{ 33 }$$

where He_n(x) are the Hermite polynomials and t_E(E, J) = E²/(Γ₂(E, J)log Λ) is the energy diffusion time. Therefore, for t > t_E/log Λ, the higher moments can be neglected, as expected from the central limit theorem.

On short timescales, the deviation from the FP propagator can be significant, as seen in Figure 2. Figure 3 shows the evolution of the energy propagator. On short timescales (t ≲ 0.01t_E), the FP approximation overestimates the contribution of small energy changes and underestimates the contribution of large energy changes to the accumulated change. Thus, formally rare events (i.e., >3σ for a Gaussian) have a much larger probability than expected by the FP treatment. For example, the probability of a 5σ Gaussian event (2.9 × 10⁻⁷) can be actually as high as 1.6 × 10⁻².

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It should be emphasized that the non-Gaussian, super-diffusive short-timescale behavior occurs in spite of the fact that the physical energy transition probability (in contrast to its simplified power-law approximation) is in fact bounded by ΔE_max, and therefore does have non-diverging moments. This seeming contradiction is reconciled by noting that on short timescales, this cutoff is irrelevant, since it is too far away in the wings of the distribution—the system does not evolve enough to "be aware" of its existence. This is a typical behavior for distributions whose variance would diverge were it not for physical cutoffs, for example, truncated Levy flights (e.g., Zumofen & Klafter 1994). Only on long enough timescales, when repeated scattering populates the wings all the way to the cutoffs, does the variance converge, the conditions for the application of the central limit theorem are satisfied, and the propagator converges to its asymptotic $\sigma \propto \sqrt{t}$ behavior, independently of the functional form of K_{E, J}(ΔE).

In Appendix B, we derive the effective untruncated (Δ_max → ∞) propagator W^s_{E, J}(Δ_tE, t). Since it includes arbitrarily strong energy exchanges, we refer to it at the "strong propagator" (SP). This propagator is given by a heavy-tailed probability distribution

$$\begin{eqnarray} W_{E,J}^{s}(\Delta _{t}E,t) & = & W_{E,J}^{0}(\Delta _{t}E,t)+\frac{W_{E,J}^{1}(\Delta _{t}E,t)}{4\log L\left(q_{2b}t\right)}, \end{eqnarray} \tag{ 34 }$$

whose leading term is a Gaussian

$$\begin{equation} W_{E,J}^{0}(\Delta _{t}E,t)=\frac{1}{\sqrt{2\pi t\Gamma _{2}^{L}(t)}}e^{-\left(\Delta _{t}E-t\Gamma _{1}^{L}(t)\right)^{2}/2t\Gamma _{2}^{L}(t)}, \end{equation} \tag{ 35 }$$

where

$$\begin{equation} \Gamma _{1,2}^{L}(t)=\Gamma _{1,2}(E,J)\log L\left(q_{2b}(E,J)t\right)\, \end{equation} \tag{ 36 }$$

and $L(x)\approx \sqrt{e^{3-2\gamma _{E}}x}$, where γ_E is the Euler constant. On long timescales, ≳ t_E/6, log L(q_2bt_E) coincides with log Λ (Equation (B15)). The second, correction term in Equation (34) is the higher order (in 1/log L(q_2bt)) even function that contributes an asymptotic tail of ∼|Δ_tE|⁻³ as |Δ_tE| → ∞,

$$\begin{eqnarray} W_{E,J}^{1}(\Delta _{t}E,t) & = & \frac{W_{c}^{1}\big(\frac{\Delta _{t}E-t\Gamma _{1}^{L}(t)}{\sqrt{t\Gamma _{2}^{L}(t)}}\big)}{\sqrt{t\Gamma _{2}^{L}(t)}}, \end{eqnarray} \tag{ 37 }$$

where W¹_c(x) is defined in Equation (B9). The next order odd correction function corresponds to the drift term (see Equation (B14)), which is small on short timescales and is omitted here. Note that unlike a Gaussian, which scales self-similarly as σ(t), the propagator in Equation (34) is not self-similar and thus does not scale globally with time. The first, central term (W⁰_{E, J}), which expresses the cumulative effect of multiple scatterings, scales with time like $\sqrt{\vphantom{A^A}\smash{\hbox{${t\Gamma _{2}^{L}(t)}$}}}\propto \sqrt{\vphantom{A^A}\smash{\hbox{${t\log (t)}$}}}$; that is, faster than $\sqrt{\vphantom{A^A}\smash{\hbox{${t}$}}}$, in a super-diffusive manner (but only logarithmically so). The second, tail term (W¹_{E, J}) expresses the effects of single rare scattering events.

On longer timescales, when the cutoffs become important and the limit Λ → ∞ is not valid, Equation (34) is no longer a good approximation of the propagator. Note that on all timescales this approximation is valid only for some central part of the propagator and it is never a good approximation of the very far tails. Therefore, this approximated propagator cannot be used in the MC simulation even if the time steps of the MC are very short.

The relaxed, steady-state configuration of a system, if such is attained on cosmologically relevant timescales, has a special significance since it is independent of the initial dynamical conditions, which are usually unknown, and since it is typically the most likely configuration to be observed. It is therefore important to understand how, and on what timescale, a system reaches steady state, and how short-timescale evolution by anomalous diffusion asymptotically approaches normal diffusion, as described by the FP approximation.

On short timescales, rare strong encounters (Δ > 1) can be neglected. However, as time progresses, even these rare events contribute significantly, and the kinematic limit (Section 3.1) must be taken into account. A key issue is to estimate their effect on the relaxation rate. In principle, the exact propagator (EP) provides a full description of all encounters. However, from a practical point of view it is difficult to implement it in an analytical or numerical scheme, due to the required triple integration (but see Goodman 1983). Instead, as we show below, the EP can be bracketed by two limiting approximations. The first, the SP W^s_E(Δ_tE) (Equation (34)), overestimates the contribution of strong encounters by taking the weak-limit transition probability (Equation (24)), and integrating it to Δ_max → ∞, i.e., beyond the weak limit (Equation (B2)). The second, the weak propagator (WP) W_E(Δ_tE), integrates only up to Δ_max = 1, where the weak limit is valid (Equation (30)). As argued below (see Figures 4 and 5), the SP, WP, and FP propagator all lead to evolution to the same steady state, on similar timescales. This then proves that so must the EP, and that the contribution of strong encounters to reaching steady state is not crucial. Table 1 summarizes the properties of the four isotropic propagators that enter our analysis of relaxation.

The inclusion of strong encounters in the transition probability can only accelerate relaxation. Therefore, the SP is expected to lead to the fastest relaxation, followed by the EP, while the WP and the FP are slowest, since they assume the weak limit. The question of the existence and values of the moments of the propagators reflect the form assumed for the transition probability, and in turn are related to the nature of the long-term evolution of the system. The weak-limit assumption introduces a cutoff on the energy exchange, which ensures that all moments exist. In particular, the neglect of all n > 2 moments assumed for the FP propagator implies that it must be a Gaussian. Both the WP and the exact one are guaranteed by the Central Limit Theorem to converge to a Gaussian. The WP by construction (Δ_max = 1) converges to the FP Gaussian. This is not the case for the EP, which takes into account the actual kinematic limits. These can be shown to correspond in the asymptotic (Gaussian) long-timescale limit to an effective Δ^eff_max ≫ 1. However, as we demonstrate below (Equation (39)), this limit is approached only on timescales much longer than the time for the stellar system to relax to steady state, and therefore has little effect on the actual time to reach steady state. This even holds for the SP (Δ_max → ∞), where all the moments diverge and the propagator never converges to a Gaussian with width ${\propto} \sqrt{t}$.

The long-term evolution of the system can be described by repeated application of the propagator on short-enough timescales. Computationally, this can be realized by an MC simulation, which follows the evolution of the orbital energy of individual stars as it changes by small random increments, drawn from the DF described by the propagator. Figure 4 shows snapshots from a set of such simulations. In each, a different propagator is used to follow the evolution of an initial δ-function distribution of test stars around an MBH, assuming a relaxed stellar background, and a reflecting boundary condition at the outer edge of the system. In all cases, the perturbation relaxes to the steady-state BW solution after t_r ∼ 10t_E. This provides numeric validation of the analytical treatment of this process (Section 2.2). As expected, the evolving DFs resulting from the WP and SP are significantly broader at early times than the FP one. This reflects their heavy-tailed probability distributions (e.g., Figure 2).

Zoom In Zoom Out Reset image size

Another demonstration of the near equivalence of the different propagators is shown in Figure 5. where the different propagators are used to follow the evolution of an initial δ-function distribution, assuming an absorbing boundary condition at the outer edge of the system. In this configuration, the time to relaxation can be defined as the time it takes the system to evaporate some large fraction of its stars. Since we are interested here mainly in comparing the evolution due to the different propagators, our conclusions do not depend on the exact fraction. For definitiveness, we choose 0.9. Both the FP and the WP drive evaporation on almost the same timescale, while the SP drives it only slightly faster, by a factor of ∼1.1 (for Δ_min = Q⁻¹ = 10⁻⁶).

Zoom In Zoom Out Reset image size

This modest acceleration in the evaporation rate by the SP, despite the inclusion of arbitrarily strong encounters, can be explained by the fact that the system evolves to its steady state much faster than the rate at which Δ > 1 events contribute to relaxation. This can be shown in a quantitative way by defining an effective, time-dependent diffusion timescale for the SP.

On long timescales, the SP is dominated by its central Gaussian, and therefore its width, σ² = tΓ₂log [L(q_2bt)] (Equation (34)), can be used to estimate an instantaneous diffusion timescale, t^ins_E(t) = E_min²/Γ₂log L(q_2bt). Therefore, at time t, the evolution of the system up to this time can be described by the FP equation with an effective (time-averaged) diffusion time

$$\begin{equation} \frac{1}{t_{E}^{\mathrm{eff}}\left(t\right)}=\frac{1}{t}\int _{0}^{t}\frac{1}{t_{E}^{\mathrm{ins}}(t^{\prime })}dt^{\prime }. \end{equation} \tag{ 38 }$$

The analytic solution of the FP equation (see Section 2.2), t_r ≈ 10t_E(E_min), suggests an association of the time for relaxation to steady state by the SP with the solution of the ansatz t_r ≈ 10t^eff_E(t_r),

$$\begin{equation} t_{r}^{{\rm SP}}\approx 10E_{\min }^{2}/[\Gamma _{2}(E_{\min })\log (\sqrt{10}e^{3/2-\gamma _{E}}/\Delta \epsilon _{\min })], \end{equation} \tag{ 39 }$$

which corresponds to an effective $\Delta \epsilon _{\max }=\sqrt{10}e^{3/2-\gamma _{E}}\approx 8$. This implies that by time t_r, the SP can be represented by an FP propagator with a Coulomb logarithm that is formally increased to ≈log (8/Δ_min) (Figure 5), thereby providing a lower limit on the time to relaxation due to the inclusion of strong encounters, t_r ≳ t^FP_rlog (1/Δ_min)/log (8/Δ_min).

This lower limit is to be compared to the time for relaxation that can be formally defined for the EP through its Γ^Λ₂ DC, where for the p = 1/4 cusp, Λ ≈ 56/Δ_min (see Equation (A16)). If the EP were to converges to the diffusive limit before steady state is attained, this would imply that t^EP_r/t_r^SP ≈ log (8/Δ_min)/log (56/Δ_min) < 1. However, this is in contradiction with t^SP_r being the lower limit on the time to relaxation, which means that the EP does not reach the diffusive limit by time t_r, and therefore cannot be fully described by an FP equation. Nevertheless, the FP equation remains a reasonable approximation for describing relaxation around an MBH, since the discrepancy in timescales is not very large. For example, for Δ_min = Q⁻¹ (Section 3.1), 0.76 < t^FP_r/t_r^SP < 0.92 in the range Q = 10³–10¹⁰ M_☉.

The propagator W_{E, J}(Δ_tE, t) (Equation (29)) is determined by the shape of the stellar cusp through the known functions Γ₁(E, J) and Γ₂(E, J) (Equations (A6) and (A7)), and by the parameters Δ_min and Δ_max. However, on the short timescales accessible through our N-body simulations, the propagator is fully determined by Γ₁, Γ₂, and the parameter Δ_min only (see Equation (34)). In practice, Γ₁(E, J) has a vanishing contribution to energy relaxation on short timescales, so that the propagator can be approximated as a functional of Γ₂ with one free parameter, Δ_min.

As shown below, careful analysis of the N-body results reveals that the propagator model provides an excellent description of the data (Figures 7 and 8). We thereby validate the functional form of Γ₂, and determine the value of Δ_min. In principle, if it were also possible to validate Γ₁ and determine Δ_max, then the full shape and range of the transition probability (Equation (24)) could be determined, and the DCs Γ^Λ_{1, 2} = Γ_{1, 2}(E, J)log Λ could be calculated reliably. As it is, we have to rely on the fact that both Γ₁ and Γ₂ were derived in the same theoretical framework, and therefore assume that the excellent fit for Γ₂ also implies a validation of Γ₁. Furthermore, since in the small deflection angle limit Δ_max ∼ O(1), setting Δ_max = 1 is expected to have only a small effect on the numerical value of log Λ = log (Δ_max/Δ_min). We find that with this empirically motivated calibration, the numerical values of the DCs lie within ≲ 20% of their theoretically predicted values.

To test and calibrate the analytic formulation of energy relaxation, we carried out a suite of Newtonian N-body simulations. Each simulation consisted of 10⁴ stars of equal mass M_⋆, and a central massive object of mass M_• = 10⁶ M_⋆. The initial semimajor axes of the stars were drawn from a ρ(a)da∝a^{2 − γ}da distribution for γ = 1.0, 1.25, 1.5, and 1.75, with random phases and isotropic orientations. The initial eccentricities were drawn from an isotropic ρ(e)de = 2ede distribution. In the Keplerian limit, these models approximately correspond to a power-law density cusp, n(r)∝r^−γ. The simulations were carried out with the parallel Nbody6++ code (implementing Hermite integration with hierarchical block time steps, the Ahmad–Cohen neighbor scheme for efficient calculation of accelerations, and the Kustaanheimo–Stiefel regularization method for dealing with close few-body sub-systems; Spurzem et al. 2001; Khalisi & Spurzem 2006).

Measuring energy relaxation in short-timescale N-body simulations is challenging, since the distribution of energy jumps is heavy-tailed, but the large energy exchange cutoff is not reached on these timescales (Section 3.2). Consequently, commonly used measures of evolution, such as the rms of the energy change, are strongly affected by rare strong exchanges and do not offer a robust characterization of the evolving energy distribution. After some experimentation, we adopted a statistically robust method, described below, whose key idea is to estimate the width of the central Gaussian component of the energy distribution as the limit of fractional moments.

Snapshots of the N-body simulation results are recorded at equally spaced times {t_i = iδt}_{i = 0, 1, 2, ...}. For each star n of the total N, the orbital energy at t_i, Eⁿ_i, is calculated and recorded. The stars are assigned to logarithmically spaced energy bins by their initial orbital energy. For each star, and for each time lag {Δt_k = kδt}_{k = 1, 2, 3, ...}, the relative energy changes are collected on non-overlapping, consecutive intervals: {Δⁿ_k} = {(E_{i + k}ⁿ − Eⁿ_i)/E_iⁿ}_{i = 0, k, 2k, 3k, ...}, thereby avoiding the complications of statistical inter-dependencies and correlations. The measured relative energy differences, {Δⁿ_k}, are drawn from the heavy-tailed parent distribution $W_{E_{i},J_{i}}(\Delta \epsilon _{k}^{n},\Delta t_{k})$, whose characteristic parameters need to be estimated from the numerical data. We use the fact that the central part of W(Δ, t) is Gaussian to a good approximation with a scale parameter σ²_ΔE(t) = tΓ₂^L(t) (Equation (34)). Since most of the data lie within the central Gaussian (Figure 8), its width, or scale parameter, can be robustly measured (Appendix B) by fractional moments of order δ < 1 (e.g., Zumofen & Klafter 1994), which give a larger weight to the central part of the distribution and a smaller one to the diverging tails,

$$\begin{eqnarray} \sigma _{\Delta E}^{2} & (t)= & \lim _{\delta \rightarrow 0}\frac{\langle |E(0)-E(t)|^{\delta }\rangle ^{2/\delta }}{2\pi ^{-1/\delta }\Gamma \left(\frac{1+\delta }{2}\right)^{2/\delta }}, \end{eqnarray} \tag{ 40 }$$

where Γ(x) is the Gamma function and we assumed that $W_{E_{i},J_{i}}(\Delta \epsilon ,t)$ is symmetric since in our case, tΓ²₁ ≪ Γ₂ (negligible drift). In particular, this relation also holds for a pure Gaussian (for any δ). We verified that the results converge rapidly as δ → 0; the results below are estimated with δ = 10⁻³.

Over the short timescales of our simulations, E and J for each star can be approximated as fixed at their initial values E_⋆ and J_⋆ for the purpose of approximating relevant parameters. Thus, using Equation (40), we define for each star a dimensionless scale parameter, which we measure directly from the data,

$$\begin{equation} \sigma _{\Delta \epsilon ,\star }^{2}=\frac{\Gamma _{2}^{L}(E_{\star },J_{\star },\tau P_{\star })\tau P_{\star }}{E_{\star }^{2}}=\lim _{\delta \rightarrow 0}\frac{\langle |\Delta \epsilon _{i,k}|^{\delta }\rangle ^{2/\delta }}{2\pi ^{-1/\delta }\Gamma \left(\frac{1+\delta }{2}\right)^{2/\delta }}, \end{equation} \tag{ 41 }$$

where Δt_k = τP_⋆, P_⋆ = P(E_⋆) and we use the minimal lag τ = 1 in order to obtain maximal statistics (Figure 6).

Zoom In Zoom Out Reset image size

We choose the energy bins to be wide enough so that the bin average over J is close enough to the isotropic average. Using Equation (41), we define the bin-averaged dimensionless scale parameter

$$\begin{eqnarray} \sigma _{\Delta \epsilon ,\mathrm{bin}}^{2}(\tau) & \equiv & \left\langle \frac{\Gamma _{2}^{L}(E_{\star },J_{\star },\tau P_{\star })\tau P_{\star }}{E_{\star }^{2}}\right\rangle _{_{{\rm bin}}}\nonumber \\ & \approx & \frac{\Gamma _{2}(E_{{\rm bin}})\tau P_{\mathrm{bin}}}{E_{\mathrm{bin}}^{2}}\log L\left(\frac{\tau P_{\mathrm{bin}}\Gamma _{2}(E_{\mathrm{bin}})}{2\Delta \epsilon _{\min }^{2}E_{\mathrm{bin}}^{2}}\right), \end{eqnarray} \tag{ 42 }$$

where E_bin and P_bin are the bin-averaged values. The second approximate equality relates the measured quantity to the theoretical model (Equation (36)), under the assumption that the weak dependence of log L(x) on x allows 〈xlog L(x)〉_bin to be approximated by 〈x〉_binlog L(〈x〉_bin). Finally, we calculate Γ₂ from Equation (D2), and fit σ²_{Δ, bin} to the N-body data to obtain the best-fit estimate of Δ_min and a goodness-of-fit score for the model (Figure 7).

Zoom In Zoom Out Reset image size

Our estimated best-fit values for Δ_min, averaged over ∼10 simulations each for various values of the power-law cusp, are presented in Table 2. We also list the relative errors resulting from the simple approximation Λ = Q, assuming Δ_max = 1. Figure 6 shows the growth of the width of the central Gaussian, σ²_{Δ, bin} as function of the time lag τ. Pure diffusion would lead to a ∝τ behavior at τ > 1. However, the results show a clear trend of a steeper than linear growth ∼τlog τ, a signature of anomalous diffusion.

Table 2. Measured Δ_min

γ	$\bar{\Delta \epsilon }_{\min }/Q^{-1}$a	$\bar{\Lambda }/Q$a, b	$\overline{\log \Lambda }/\log Q$ a,b	$\overline{\chi _{\mathrm{red}}^{2}}(\bar{\Delta \epsilon }_{\min })$	nsim
1.0	5.3 ± 0.5	0.19 ± 0.02	0.88 ± 0.01	1.2 ± 0.3	8
1.25	6.7 ± 0.7	0.15 ± 0.02	0.86 ± 0.01	1.3 ± 0.2	12
1.5	9.2 ± 1.3	0.11 ± 0.02	0.84 ± 0.01	1.5 ± 0.5	9
1.75	13.7 ± 3.6	0.07 ± 0.02	0.81 ± 0.02	2.0 ± 0.6	9

Notes. ^aFor Q = 10⁶, best-fit values for N-body results (Section 4.2). ^b$\bar{\Lambda }=\Delta \epsilon _{\max }/\bar{\Delta \epsilon }_{\min }$, where Δ_max = 1.

Download table as:  ASCIITypeset image

Finally, to verify our analysis we plot the histogram of the orbital energy jumps Δ_{i, k} in a given energy bin and compare it to the bin-averaged propagator

$$\begin{eqnarray} W_{E_{\mathrm{bin}}}(\Delta \epsilon ,t) & \simeq & \langle W_{E_{\star },J_{\star }}(\Delta \epsilon ,t)\rangle _{E_{\mathrm{bin}}}, \end{eqnarray} \tag{ 43 }$$

where we evaluated $W_{E_{\star },J_{\star }}(\Delta \epsilon ,t)$ numerically for each star using the fitted ΔE_min and Equation (A7). An example is shown in Figure 8. For comparison we also show the expected bin-averaged FP propagator $\langle W_{E_{\star },J_{\star }}^{{\rm FP}}(\Delta \epsilon ,t)\rangle$ and the bin-averaged central Gaussian of the short-timescale propagator (Equation (34)).

Zoom In Zoom Out Reset image size