THE EFFICIENCY OF RESONANT RELAXATION AROUND A MASSIVE BLACK HOLE

The NR time for E relaxation, T^E_NR, corresponds to the time it takes noncoherent two-body interactions to change the stellar orbital energy by order of itself, |ΔE| ∼ E (by stellar dynamical definition convention, E>0 for a bound orbit). Similarly, the NR time for J relaxation, T^J_NR, corresponds to the time it takes the stellar orbital angular momentum to change by order of the circular angular momentum |ΔJ| ∼ J_c, where near an MBH of mass M, $J_{c} = GM/\sqrt{2E}$. The E-relaxation timescale can be estimated by considering the rate Γ of gravitational collisions at a relative velocity v between a test star and N field stars of mass m with a space density n ∼ N/R³ in a system of size R, at the minimal impact parameter where the small angle deflection assumption still holds, r_min ∼ Gm/v². The collision rate is then Γ ∼ nvr²_min ∼ G²m²n/v³. Taking into account also collisions at larger impact parameters increases the rate by the Coulomb logarithm factor ln Λ ∼ ln(R/r_min), so that T^E_NR ∼ v³/(G²m²nln Λ). Near the MBH v² ∼ GM/R, and so ln  Λ ∼ ln  Q, where Q ≡ M/m is the mass ratio. When the stars move under the influence of the central MBH (Q ≫ N), the relaxation time can be expressed as T^E_NR ∼ (M/m)²P/(Nln Λ), where $P = 2\pi \sqrt{R^{3}/GM}$ is the Keplerian period.

Rauch & Tremaine (1996, RT96) parameterized the uncertainties in the evolution of the orbital E, J, and J by a set of numeric coefficients (α_Λ, η_s,v, β_s,v), whose operational definition is tied to the procedure by which they are estimated (left unspecified by RT96). Here, we adopt this notation and estimate these coefficients by measuring the rms change in energy and angular momentum over the stellar population at binned time lags (Section 3.2). The NR changes in E, J, and J over the dimensionless time lag τ ≡ (t₂ − t₁)/P₁ are thus formally defined as²

$$\begin{eqnarray} \delta E(\tau) & \equiv & \langle (|E_{2} - E_{1}|/E_{1})^{2}\rangle ^{1/2} = \alpha _{\Lambda }\sqrt{N}(m/M)\sqrt{\tau }\,, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \delta J(\tau) & \equiv & \langle (|J_{2} - J_{1}|/J_{c,1})^{2}\rangle ^{1/2} = \eta _{s\Lambda }\sqrt{N}(m/M)\sqrt{\tau }\,, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \delta \mathbf {J}(\tau) & \equiv & \langle (|\mathbf {J}_{2} - \mathbf {J}_{1}|/J_{c,1})^{2}\rangle ^{1/2} = \eta _{v\Lambda }\sqrt{N}(m/M)\sqrt{\tau }, \end{eqnarray} \tag{ 3 }$$

where 〈 ⋅ ⋅ ⋅ 〉 designates the average over all stars in the time-lag bin τ, $\alpha _{\Lambda } \equiv \alpha \sqrt{\hbox{ln}\,\Lambda }$ and $\eta _{s,v\Lambda } \equiv \eta _{s,v}\sqrt{\hbox{ln}\,\Lambda }$ are dimensionless constants, whose exact values are system dependent. Accurate determination of their values requires detailed calculations or simulations. The corresponding NR timescales are related to these coefficients by T^E_NR = (M/m)²P/(Nα²_Λ), T^J_NR = (M/m)²P/(Nη²_sΛ), and T^J_NR = (M/m)²P/(Nη²_vΛ).

When the gravitational potential has symmetries that restrict the orbital evolution, for example, to fixed ellipses in the potential of a point mass, or to planar annuli in a spherical potential, the perturbations on a test star are no longer random, but correlated. This leads to a coherent changes ΔJ = Tt on times P ≪ t < t_ω by the residual torque $\left|\mathbf {T}\right| \sim \sqrt{N}Gm/R$ exerted by the N randomly oriented, orbit-averaged mass distributions of the surrounding stars (mass wires for elliptical orbits in a Kepler potential, mass annuli for rosette-like orbits in a spherical potential). The coherence time t_ω is set by deviations from perfect symmetry, which lead to a gradual orbital drift and to the randomization of T. For example, the enclosed stellar mass leads to non-Keplerian retrograde precession; general relativity leads to prograde precession. Ultimately, the coherent torques themselves randomize the orbits. The effective coherence time is set by the shortest decoherence (quenching) process. The accumulated change over t_ω, |ΔJ_ω| ∼ |Tt_ω|, then becomes the basic step-size, or mean free path in J space, for the long-term (t ≫ t_ω) random-walk phase ($\propto \sqrt{t}$) of RR of J. Since this step-size is large, RR can be much faster than NR. The RR timescale T_RR is then defined by $\left|\Delta J\right|/J_{c} = (\left|\Delta J_{\omega }\right|/J_{c})\sqrt{t/t_{\omega }} \equiv \sqrt{t/T_{{\rm RR}}}$. Note that the relaxation of E is not affected by RR because the potential of the system is stationary on the coherence timescale, and so E changes incoherently on all timescales. The torques exerted by elliptical mass wires in a Kepler potential can change both the direction and magnitude of J. In contrast, the torques exerted by planar annuli can only change the direction of J. The long-term random-walk $\sqrt{\tau }$ growth continues until |ΔJ|/J_c ∼ O(1), at which point |ΔJ|/J_c and |ΔJ|/J_c can no longer grow, but J does continue to change its value randomly.

Here, we consider only Newtonian dynamics. The coherence timescale for scalar RR is determined by the time it takes for the orbital apsis to precess by angle ∼π due to the potential of the enclosed stellar mass ("mass precession"),

$$\begin{equation} t_{\omega }=t_{M}=A_{M}\left(M/Nm\right)P\,, \end{equation} \tag{ 4 }$$

where A_M is an O(1) factor reflecting the approximations in this estimate. The coherence timescale for vector RR is determined by the time it takes for the coherent torques to change by |ΔJ| ∼ J_c (self-quenching³),

$$\begin{equation} t_{\omega }=t_{\phi }=A_{\phi }(\sqrt{N}/\mu)P\simeq \big[A_{\phi }(M/m)\big/\sqrt{N}\big]P\,, \end{equation} \tag{ 5 }$$

where A_ϕ is an O(1) factor,⁴ μ = Nm/(M + Nm) and where the approximate equality is for the Keplerian limit Nm ≪ M. The RR changes in J and J during the coherent phase (τ < τ_ω) are defined in terms of the rms change as

$$\begin{eqnarray} \delta J(\tau) & \equiv & \langle (|J_{2}-J_{1}|/J_{c,1})^{2}\rangle ^{1/2}=\beta _{s}\sqrt{N}(m/M)\tau, \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} \delta \mathbf {J}(\tau) & \equiv & \langle (|\mathbf {J}_{2}-\mathbf {J}_{1}|/J_{c,1})^{2}\rangle ^{1/2}=\beta _{v}\sqrt{N}(m/M)\tau, \end{eqnarray} \tag{ 7 }$$

where the dimensionless coefficients β_s and β_v depend on the parameters of the system and reflect the uncertainties introduced by the various approximations and simplification of this analysis. Accurate determination of their values requires detailed calculations or simulations.

The scalar RR change in J on time lags τ ≫ τ_M is then

$$\begin{equation} \delta J(\tau) \equiv \langle (|J_{2} - J_{1}|/J_{c,1})^{2}\rangle ^{1/2} = \beta _{s}\sqrt{A_{M}(m/M)}\sqrt{\tau }, \end{equation} \tag{ 8 }$$

and the scalar RR timescale is

$$\begin{equation} T_{{\rm RR}}^{J}=\big[(M/m)/A_{M}\beta _{s}^{2}\big]P. \end{equation} \tag{ 9 }$$

The RR efficiency factor χ = (β_s/β_s,RT96)² defined by Hopman & Alexander (2006) for an assumed value A_M = 1 (Section 3.3), expresses how much shorter T^J_RR is relative to the value estimated by RT96. Scalar RR is faster than NR by a factor T^J_NR/T^J_RR ∝ (M/m)/Nln  Λ. Similarly, the vector RR change in J on time lags τ ≫ τ_ϕ in the Keplerian limit is

$$\begin{equation} \delta \mathbf {J}(\tau) \equiv \langle (|\mathbf {J}_{2} - \mathbf {J}_{1}|/J_{c,1})^{2}\rangle ^{1/2} = \beta _{v}\sqrt{A_{\phi }\sqrt{N}(m/M)}\sqrt{\tau }, \end{equation} \tag{ 10 }$$

and the vector RR timescale is

$$\begin{equation} T_{{\rm RR}}^{\mathbf {J}}=\big[(M/m)/\sqrt{N}A_{\phi }\beta _{v}^{2}\big]P. \end{equation} \tag{ 11 }$$

RT96 performed a limited set of near-Keplerian simulations to check their predictions. They analyzed the results both star by star and in the average and observed the coherent growth of |ΔJ|/J_c and |ΔJ|/J_c relative to the simulation's initial values. Although the evolution of these quantities for any single star was very noisy and the proportionality factors had a very large scatter, RR was clearly observed, as predicted.

Our N-body code uses a 5th order Runge–Kutta integrator with individual time steps and optional pairwise KS regularization (Kustaanheimo & Stiefel 1965), without gravity softening. The time steps were chosen to conserve total energy at the level of |ΔE_tot|/E_tot ∼ O(10⁻⁵)–O(10⁻⁸), well below the NR energy changes expected in the simulations.

To simulate RR reliably, the N-body code must maintain small-enough numerical deviations from the conserved Keplerian symmetries. These could lead to numerical quenching of RR by a drift of the orbital apsis or of the direction of the angular momentum, which change the torque felt or exerted by the orbit, and to a lesser degree by a drift in the orbital period (equivalently, the energy), which changes the shape of the orbit, and hence its moment of inertia. We estimated these numerical drifts by integrating, over many orbits, a highly eccentric two-body system (e = 0.99995, Q = 3 × 10⁶), since most of the numerical errors accumulate at periapse, where the acceleration is largest. We tracked the change in the period, ΔP = (P − P₀)/P₀, in the angular momentum ΔJ = |J − J₀|/J₀ and in the direction of the apsis $\Delta \mathbf {\hat{e}} = |\mathbf {\hat{e}}-\mathbf {\hat{e}}_{0}|/|\mathbf {\hat{e}}_{0}|$, where $\hat{\mathbf {e}} = \mathbf {e}/|\mathbf {e}|$ and e = v ⊗ (r ⊗ v)/(M + m) − r/|r| is the Laplace–Runge–Lenz vector. We find that after τ_sim = 4.5 × 10⁴, ΔP ∼ 8 × 10⁻⁵, ΔJ ∼ 3 × 10⁻⁴, and $\Delta \mathbf {\hat{e}} \sim 9\times 10^{-8}$, which correspond to drift timescales of τ_ΔP ∼ (π/ΔP)τ_sim = 2 × 10⁹, τ_ΔJ ∼ 5 × 10⁸, and $\tau _{\Delta \mathbf {\hat{e}}} \sim 3\times 10^{12}$. These are many orders of magnitude longer than our longest simulations (e.g., Figure 5), and can therefore be safely neglected.

We carried out two sets of simulations. The first, designed to study the short-term behavior of RR (τ_sim < τ_M) and to measure NR and RR parameters α, η_s,v, and β_s,v, consisted of 200 particles (including the MBH as a fixed particle, since testing showed that allowing the MBH to move did not change the results significantly, but slowed down the simulations substantially by requiring smaller time steps and KS regularization). The initial orbital semimajor axes were randomly drawn from a ρ(a)da ∝ a^2−γda distribution for γ = 1,  1.5,  and 1.75, for a in the range (0,  1/2), with eccentricities drawn from a thermal ρ(e)de = 2ede distribution, with random phases, orbital orientations, and isotropic velocities. This distribution approximates an r^−γ number density distribution with an outer cutoff at radius r = R = 1 from the MBH (in dimensionless units where G = 1, M + Nm = 1). These stellar cusps span a wide range of possible physical scenarios (e.g., Bahcall & Wolf 1977), and in particular those of LISA targets, which are expected to be relaxed galactic nuclei (Alexander 2007). A typical simulation lasted a few × 100 system orbital times and resulted in few × 100 snapshots of the system configuration. In order to decrease the statistical errors, we ran n_sim = 12–15 simulations with different initial conditions for each of the (γ, Q) models we studied.

The second set of simulations was designed to study the long-term behavior of RR (τ_sim ≫ τ_ϕ), the transition from the coherent RR phase to the random-walk phase, and the asymptotic saturation limits. Due to the high computational load, this set (n_sim = 7) was limited to a relatively small mass ratios, Q ⩽ 10⁴, smaller number of particles (50 ⩽ N ⩽ 200), and flat cusps with γ = 1, where time-consuming close star-star interactions are rarer. A typical simulations lasted few × 10⁴ system orbital times (τ_sim ∼ 10⁴).

Given the high computational cost of the N-body simulations, and the very large variance in the evolution of individual orbits, it is essential to extract the RR signal as efficiently and robustly as possible from the simulated data. After some experimentation, we adopted the correlation analysis for detecting and measuring RR in N-body simulation snapshots. Our procedure is as follows.

• 1.  
  The stellar phase space coordinates are transformed to the rest frame of the MBH, which is almost identical with the center of mass in the near-Keplerian system close to the MBH.
• 2.  
  The energy (E⁽ⁿ⁾_i), angular momentum (J⁽ⁿ⁾_i), circular angular momentum (J⁽ⁿ⁾_c,i), and Keplerian period (P⁽ⁿ⁾_i) are calculated for the nth star (n = 1...N) at discrete times t_i in the simulation ("snapshots").
• 3.  
  A normalized time lag τ⁽ⁿ⁾_ji = (t_j − t_i)/P⁽ⁿ⁾_i is assigned to each pair of times (t_j>t_i). For each of the N stars, we calculate the normalized energy and angular momentum differences at all lags, (|ΔE|/E)⁽ⁿ⁾_ji = |E⁽ⁿ⁾_j − E⁽ⁿ⁾_i|/E⁽ⁿ⁾_i, (|ΔJ|/J_c)⁽ⁿ⁾_ji = |J⁽ⁿ⁾_j − J⁽ⁿ⁾_i|/J⁽ⁿ⁾_c,i, and (|ΔJ|/J_c)⁽ⁿ⁾_ji = |J⁽ⁿ⁾_j − J⁽ⁿ⁾_i|/J⁽ⁿ⁾_c,i.
• 4.  
  The differences from all stars are binned into discrete τ bins according to their associated τ_ij. The bin rms values and their standard deviations are calculated and plotted against the bin's average time lag τ = 〈τ_ij〉, thereby creating the correlation curve.

By using all possible time lags recorded in the simulation, this approach makes maximal use of the entire data set and averages over the strongly fluctuating individual relaxation curves (e.g., the RT96 procedure, Section 2.2). However, this method is not entirely free of bin-to-bin bias. Since the number of orbital periods completed by a star with a mean period P⁽ⁿ⁾_i over the simulation time t_sim is t_sim/P⁽ⁿ⁾ = τ⁽ⁿ⁾_max, long-period stars will not contribute to a high-τ bin, 〈τ〉_i if τ⁽ⁿ⁾_max < 〈τ〉_i. Conversely, short-period stars will not contribute to a low-τ bin, 〈τ〉_i, if τ⁽ⁿ⁾_min>〈τ〉_i, where τ⁽ⁿ⁾_min = min_j(t_j+1 − t_j)/P⁽ⁿ⁾ is set by the minimal time difference between consecutive snapshots. In extreme cases, some stars may not contribute to the relaxation curve at any 〈τ〉_i. Since long- and short-period stars could well have systematically different responses to RR, this introduces bias to the relaxation curve. This bias can be reduced, at the cost of losing some information, by using only the middle range of the τ bins, or at a computational cost, by longer simulations with a higher snapshot rate.

The correlation curves δJ(τ) and δJ(τ) reflect the joint effects of NR and RR. The coefficients η_s,v and β_s,v are formally defined and measured by fitting the measured correlation curves in the coherent phase to the functions

$$\begin{eqnarray} \delta J(\tau) & = & \sqrt{N}(m/M)\sqrt{\eta _{s\Lambda }^{2}\tau +\beta _{s}^{2}\tau ^{2}}\,, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \delta \mathbf {J}(\tau) & = & \sqrt{N}(m/M)\sqrt{\eta _{v\Lambda }^{2}\tau +\beta _{v}^{2}\tau ^{2}}\,, \end{eqnarray} \tag{ 13 }$$

where the two terms in the square root express the contributions of NR and RR, respectively. Note that on short timescales, τ < τ₀ ≡ (η_sΛ/β_s)², the RR correlation curve rises as $\sqrt{\tau }$ due to the effect of NR. It then rises as τ in the RR-dominated coherent phase at times τ₀ < τ < τ_ω, before turning over again to a $\sqrt{\tau }$ rise in the accelerated random-walk phase at τ>τ_ω (this last phase is not modeled by these fitting functions).

Figure 1 shows the measured correlation curves in a typical simulation for τ ≪ τ_ω: δE(τ), δJ(τ), and δJ(τ). The NR and RR coefficients were derived from the best two-parameter fits of Equations (1), (12), and (13) to the data. To control the star to star variance, we limited our analysis to time-lag bins that sampled at least 0.75 of the stars in the simulation. The excellent fit of the data points to the predicted correlation curves seen in Figure 1 indicates that RR is present and measurable.

Zoom In Zoom Out Reset image size

To quantify the transition between the coherent and random-walk scalar RR phases and measure the coefficient A_M that sets the scalar RR coherence time (Equation (4)), and therefore the scalar RR timescale (Equation (9)), we fitted the scalar RR correlation curve up to τ_ϕ by the composite function

$$\begin{equation} \delta J_{b}(\tau ;\tau _{b})=\left\lbrace \begin{array}{@{}lr} \delta J(\tau) & \tau <\tau _{b}\\ \ms \delta J(\tau _{b})\sqrt{\tau /\tau _{b}} & \tau \ge \tau _{b},\end{array}\right. \end{equation} \tag{ 14 }$$

where δJ(τ) (Equation (12)) is the best fit of the short timesscale NR and coherent RR phases up to τ_M and τ_b is the break timescale. The best-fit coefficient A_M is then given by the best-fit break time lag, A_M = (N/Q)τ_b.

Although the correlation analysis stabilizes against star to star scatter in a single simulation, we still find a large simulation to simulation scatter in the derived values of the coefficients. We therefore constructed a large grid of near-Keplerian, short-timescale (τ₀ ≪ τ_sim < τ_M) models, where coherent RR should be clearly detected.

We summarize our results in Table 1 and in Figures 2–4. The coefficients α_Λ and η_s,vΛ do not directly express the intrinsic properties of NR, since they should vary as $\sqrt{\hbox{ln}\,\Lambda } \simeq \sqrt{\ln \,Q} \simeq 4(1\,\pm\,0.07)$ over the Q range of our models. This small fractional difference and the relatively low statistics of our simulation suite may explain why we do not detect a clear Q dependence in these quantities. Since we do not see a dependence on Q, and the γ dependence seems weak (α_Λ shows hints of an increase with γ, and η and β of a decrease with γ), we adopt as the best-fit estimates of the values of α, η_s,v, and β_s,v their grand average over all the simulations (Table 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Measured NR and RR Coefficients^a

γ	Q	nsim	$\bar{\alpha }_{\Lambda }$	$\bar{\eta }_{s\Lambda }$	$\bar{\eta }_{v\Lambda }$	$\bar{\beta }_{s}$	$\bar{\beta }_{v}$
1	106	12	10.48 ± 0.43	4.58 ± 0.38	7.31 ± 0.61	1.14 ± 0.07	1.95 ± 0.08
1	107	12	9.37 ± 0.32	4.18 ± 0.26	6.69 ± 0.46	1.13 ± 0.07	1.98 ± 0.10
1	108	12	10.27 ± 0.36	4.26 ± 0.22	7.02 ± 0.47	1.13 ± 0.06	1.89 ± 0.08
1.5	106	13	12.61 ± 0.90	4.12 ± 0.48	6.65 ± 0.61	1.00 ± 0.05	1.77 ± 0.07
1.5	107	12	11.35 ± 0.46	4.09 ± 0.35	6.58 ± 0.57	1.05 ± 0.06	1.86 ± 0.07
1.5	108	12	10.51 ± 0.51	3.99 ± 0.21	6.57 ± 0.49	1.09 ± 0.07	1.92 ± 0.09
1.75	106	14	12.06 ± 0.54	3.97 ± 0.35	6.50 ± 0.48	0.95 ± 0.06	1.68 ± 0.06
1.75	107	12	12.29 ± 0.43	3.39 ± 0.14	5.29 ± 0.18	0.99 ± 0.07	1.75 ± 0.09
1.75	108	15	11.90 ± 0.43	3.89 ± 0.39	6.20 ± 0.52	1.01 ± 0.05	1.73 ± 0.05
Grand averageb			〈α〉	〈ηs〉	〈ηv〉	〈βs〉	〈βv〉
		114	2.82 ± 0.05	1.01 ± 0.03	1.64 ± 0.04	1.05 ± 0.02	1.83 ± 0.03

Notes. ^aThe quoted errors are the errors on the mean (the sample rms is $\sqrt{n_{{\rm sim}}}$ times larger). ^b$\langle \alpha \rangle = \langle \alpha _{\Lambda }/\sqrt{\ln \,\Lambda }\rangle$, $\langle \eta _{s,v}\rangle = \langle \eta _{s,v\Lambda }/\sqrt{\ln \,\Lambda }\rangle$ over all simulations for Λ = Q.

Download table as:  ASCIITypeset image

To explore the long-term behavior of the system, verify that the short-timescale coherent torquing phase leads to long-term relaxation (rather than oscillatory behavior, for example) and to estimate the coefficients that set the coherence timescales, we also conducted several very long simulations (τ ∼ 10⁴ ≫ τ_ϕ). The resulting RR correlation scalar curves display a first clear detection of a transition from the asymptotically ∝τ coherent RR phase to the ∝$\sqrt{\tau }$ RR relaxation phase (Figure 5) at the mass precession coherence time τ_M = A_MQ/N (Equation (4)). The mean best-fit estimate (Equation (14)) over seven long-timescale simulations is A_M = 0.99 ± 0.06, where the error is dominated by the simulation to simulation scatter. As discussed further below (Section 4.2), the vector RR curve does not display a coherent rise beyond τ_M, since by its definition it includes a scalar component, and neither does it display a distinct random-walk phase beyond τ_ϕ, since the transition to the asymptotic saturation limit occurs soon after, at τ ≳ τ_ϕ. This behavior makes it harder to quantify and measure A_ϕ by a fit to the correlation curve, and we do not attempt it here.

Zoom In Zoom Out Reset image size

The relations between the correlation functions of scalar and vector relaxation, parameterized by the coefficients η_s,v and β_s,v (Equations (12) and (13)), reflect the symmetries and dimensionality of the torquing process. In the limit |ΔJ|/J ≪ 1, the ratios η_v/η_s and β_v/β_s can be expressed as (〈|ΔJ²|〉/〈ΔJ²_∥〉)^1/2 = (1 + 〈ΔJ²_⊥〉/〈ΔJ²_∥〉)^1/2, where ΔJ_∥ is the change along J and ΔJ_⊥ perpendicular to it, in the orbital plane.

In the case of NR relaxation, and in the limit of close (impulsive) encounters, the torquing is two-dimensional, 〈ΔJ²_⊥〉 = 〈ΔJ²_∥〉, and $\eta _{v}/\eta _{s} = \sqrt{2}$. This can be seen by expressing dJ to leading order as dJ/J_c ≃ (dr/r) ⊗ (v/v) + (r/r) ⊗ (dv/v), where r and v are the instantaneous position and velocity of the scattered star at closest approach to the scatterer. The virial velocity in a self-gravitating system of typical size R = 1 is $V = \sqrt{GM/R} = 1$ (for G = M = 1, m ≪ M) and the crossing time is t_c ∼ R/V = 1. The acceleration in a two-body encounter with impact parameter b is a ∼ Gm/b² ∼ m/b², and the encounter time is t_b ∼ b/V ∼ b. The typical change in the velocity of the interacting stars is then dv ∼ at_b ∼ m/b, and in their position dr ∼ at²_b ∼ m. For a typical orbit with r ∼ R = 1 and v ∼ V = 1, the ratio of relative changes due to the encounter is (dv/v)/(dr/r) ∼ 1/b. Therefore, in an impulsive, local scattering event where b ≪ R = 1, dJ/J_c ∼ (r/r) ⊗ (dv/v) is two-dimensional. In practice, however, longer range nonimpulsive encounters that last a substantial fraction of the orbit also contribute to NR. In such encounters dr/r is not negligible, dr ∼ ∫dt'∫dt a(t) and dv ∼ ∫dt a(t) are not colinear, and r(t) and v(t) change in the course of interaction. Therefore, NR is not expected to be exactly two-dimensional, and we indeed find in our simulations that 〈η_v〉/〈η_s〉 = 1.61 ± 0.06 lies between $\sqrt{2}$ and $\sqrt{3}$.

By contrast, in the case of RR the torquing is global, averaged over the entire orbit and exerted over many orbital periods, and the impulsive encounter argument for two-dimensional torquing does not apply. We find in our simulations that 〈ΔJ²_⊥〉 ≃ 2〈ΔJ²_∥〉, when averaged over a thermal eccentricity distribution, and that $\left\langle \beta _{v}\right\rangle /\left\langle \beta _{s}\right\rangle =1.74\pm 0.04 \simeq \sqrt{3} = 1.73$. The difference from 〈η_v〉/〈η_s〉 is statistically significant, and presumably reflects the qualitatively different nature of NR and RR. However, the result 〈β_v〉/〈β_s〉 ≃  $\sqrt{3}$ is not due to isotropic torques. Rather, the decomposition of ΔJ along the normalized axes of the orbital ellipse ($\hat{\mathbf {e}},\,\hat{\mathbf {\mathbf {b}}},\,\mathbf {\hat{J}}$), where $\hat{\mathbf {\mathbf {b}}}$ is in the direction of the semiminor axis, reveals that 〈ΔJ²_b〉 ∼ 2〈ΔJ²_∥〉 ≫ 〈ΔJ²_e〉. The small values of ΔJ_e can be understood if the residual force F acting on the mass wire is approximately spatially constant, since then the ellipse's center of mass lies along e, and the torque $\mathbf {T} \propto \mathbf {\hat{e}}\otimes \mathbf {F}$ has no component along $\hat{\mathbf {e}}$ (e.g., Landau & Lifshitz 1976, Section 34). However, the finding 〈ΔJ²_b〉 ∼ 2〈ΔJ²_∥〉, and the resulting coincidental similarity to isotropic torquing when averaged over a thermal eccentricity distribution, remain to be explained.

The long-term behavior of the correlation curves is quite distinct from its short-term (τ < τ_M) behavior. Two regimes are apparent: the intermediate timescale regime (τ_M < τ < τ_ϕ), when J grows by random walk, but the direction of the angular momentum $\hat{\mathbf {J}}$ is still changing linearly with time, and the long timescale, when both the scalar and vector correlation curves reach their asymptotic saturation limits.

Intermediate timescale (τ_M < τ < τ_ϕ). The vector RR correlation curve does not explicitly display a coherent (∝τ) rise beyond τ_M, even though the coherence time for vector RR τ_ϕ ≫ τ_M. This results from the fact that as defined, δJ describes the change in both the direction and magnitude of the angular momentum. On the intermediate timescale, τ_M < τ < τ_ϕ, the scalar component in δJ evolves as $\sqrt{\tau }$, while the vector component evolves as τ. This results in a mixed behavior that is intermediate between coherent growth and random walk. This can be seen explicitly by decomposing (δJ)² = δJ² + 〈2(J²/J²_c)[1 − cos(Δθ)]〉, where the approximations J_c = const (valid since ΔE is small for τ_ϕ ≪ τ ≪ τ_NR over a wide range, $\tau _{{\rm NR}}/\tau _{\phi } = Q/(\sqrt{N}\alpha ^{2}\,\log \,Q) \gg 1$), and J₁ = J₂ = J are assumed, and where cos Δθ = J₁ · J₂/J². The angle Δθ is related to δJ by $J_{c}(\delta \mathbf {J})/J = \sqrt{2(1-\cos \,\Delta \theta)} = \beta _{v}(\sqrt{N}/Q)\tau$ for τ < τ_ϕ (Equation (7)), and δJ = τ/τ_sRR for τ>τ_M (Equation (9)). It then follows that $\delta \mathbf {J} = (\beta _{s}A_{M}/\sqrt{N})\sqrt{(\tau /\tau _{M})+(\beta _{v}/\beta _{s})^{2}(\tau /\tau _{M})^{2}}$ for τ>τ_M. This mixed power-law behavior further modified as τ → τ_ϕ, δJ → 1 and the correlation curve approaches its asymptotic saturation limit.

Long timescale (τ_ϕ < τ). When the time lag τ ≫ τ_ϕ ≫ τ_M, both J and e are completely randomized and the change in the angular momentum saturates at δJ(τ ≫ τ_ϕ) = 1/3 and δJ(τ ≫ τ_ϕ) = 1 (Figure 5). This can be seen by considering the population rms change in the angular momentum between times τ₁ and τ₂, δJ = 〈(|J₂ − J₁|/J_c)²〉^1/2 and δJ = 〈(|J₂ − J₁|/J_c)²〉^1/2, with the approximation J_c = const. For the isotropic thermal angular momentum distribution modeled here, the parent probability density of J is ρ_J = 2J/J²_c for J = 0 to J_c. Assuming no correlation between J₁ and J₂ after long-enough time lags, and replacing the population rms by the relevant moments J^a₁J^b₂ over ρ_J, $\overline{(J_{1}/J_{c})^{a}(J_{2}/J_{c})^{b}} = 4/[(2+a)(2+b)]$, yields $\delta J = \big(\overline{J_{1}^{2}} - 2\overline{J_{1}J_{2}} + \overline{J_{2}^{2}}\big)^{1/2}/J_{c} = 1/3$ and $\delta \mathbf {J}=\big(\overline{J_{1}^{2}} - 2\overline{\mathbf {J}_{1}\,{\bm \cdot}\, \mathbf {J}_{2}} + \overline{J_{2}^{2}}\big)^{1/2}/J_{c} = 1$, as measured.

Our measured values for β_s,v are in excellent agreement with the results of Gürkan & Hopman (2007), who modeled the torquing stars by the static torque field of 10⁴ random fixed ellipses, and investigated the dependence of the RR efficiency on the eccentricity of test orbits. They report⁵ β_s(e)/2π = 0.25e and β_v/2π = 0.28(0.5 + e²), which for the thermal eccentricity distribution n(e)de = 2e modeled here, translate to 〈β_s〉 = 1.05 and 〈β_v〉 = 1.76 (e.g., Table 1). We note that the advantage of the full N-body approach used here, beyond verifying the approximations and assumptions that go into the fixed torques approach (orbital averaging, neglect of close encounters) is that it enables the study of the transition from and to random-walk phases, which cannot be reproduced by a fixed background. Specifically, the transition from the early NR-dominated phase to the coherent RR phase, and the subsequent transition out of the coherent RR phase to the RR random-walk phase due to quenching (especially in the case of the self-quenched vector RR, Section 2.2).

RT96 derived from their simulations smaller mean values for β_s (0.53, as compared to 1.05 here) and for α_Λ and η_sΛ (3.09 and 1.37 as compared to 11.25 and 4.05 here). At least part of the difference in α_Λ and η_sΛ can be traced to their use of a softening length = 10⁻²R ≫ Gm/v² ∼ (m/M)R in the calculation of the gravitational force. This decreases the effective value of the Coulomb factor, since Λ ∼ R/max(, Gm/v²) = 10² and so $\sqrt{\ln \,\Lambda } \simeq 2.1$, about twice as small as in our simulations. Indeed, RT96 noted that a decrease in the softening length to ε = 10⁻⁴R led to an increased value of α_Λ = 5.5 ± 0.2. It is difficult to trace the specific reason for the discrepancy between our best estimate value for β_s and that derived by RT96, given the many differences in both the simulations and the methods of analysis. However, we note that the statistics here are better due to the larger number of simulations, more efficient use of the data, and the more rigorous analysis.

We conclude the discussion of our results by briefly considering the implications of this revised value of β_s for EMRI rates. Hopman & Alexander (2006) assumed A_M = 1, parameterized the RR efficiency by a factor χ = [β_s/β_s,RT96]², and derived the χ dependence of the branching ratio of the EMRI and infall (plunge) rates (Figure 6). The RT96 value χ = 1 happens to lie very close to the maximum of the RR-accelerated EMRI rate. We confirm here that A_M = 0.99 ±  0.06 and find β_s/β_s,RT96 = 1.98 ± 0.04, which corresponds to a factor 4.9 increase in the EMRI rate compared to that estimated for NR only, but is a factor ∼1.5 smaller than implied by the RT96 value, because the higher RR efficiency leads to a higher plunge rate at the expense of the inspiral rate.

Zoom In Zoom Out Reset image size