GRAVITATIONAL ENCOUNTERS AND THE EVOLUTION OF GALACTIC NUCLEI. III. ANOMALOUS RELAXATION

To account for anomalous relaxation, the angular-momentum diffusion coefficients of Equation (13) are modified as follows:

$$\begin{eqnarray}\langle {\rm{\Delta }}{\mathcal{R}}\rangle & = & \langle {\rm{\Delta }}{\mathcal{R}}{\rangle }_{\mathrm{CK}}+{g}_{1}({\mathcal{E}},{\mathcal{R}})\langle {\rm{\Delta }}{\mathcal{R}}{\rangle }_{\mathrm{RR}},\\ \langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle & = & \langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{CK}}+{g}_{2}({\mathcal{E}},{\mathcal{R}})\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{RR}}.\end{eqnarray} \tag{ 24 }$$

The functions ${g}_{1}({\mathcal{R}})$ and ${g}_{2}({\mathcal{R}})$ should have certain properties. Both g₁ and g₂ should tend to one as ${\mathcal{R}}\to 1$. Since

$$\begin{eqnarray*}&&\langle {\rm{\Delta }}{\mathcal{R}}{\rangle }_{\mathrm{RR}}\to 2A({\mathcal{E}}),\ \langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{RR}}\to 4A({\mathcal{E}}){\mathcal{R}}\end{eqnarray*}$$

as ${\mathcal{R}}\to 0$ (Equation (14)), we require ${g}_{1}\to {{\mathcal{R}}}^{2}$ and ${g}_{2}\to {{\mathcal{R}}}^{2}$ for small ${\mathcal{R}}$ so that the small-${\mathcal{R}}$ behavior of Equation (22) is reproduced. The transition between the two regimes should occur at ${\mathcal{R}}\sim {{\mathcal{R}}}_{\mathrm{SB}}$ for both functions.

An ad hoc functional form for g₂ that satisfies these requirements is

$$\begin{eqnarray}&&{g}_{2}({\mathcal{E}},{\mathcal{R}})={\{1+{[\displaystyle \frac{{{\mathcal{R}}}_{2}({\mathcal{E}})}{{\mathcal{R}}}]}^{n}\}}^{-2/n}\end{eqnarray} \tag{ 25 }$$

where ${{\mathcal{R}}}_{2}\approx {{\mathcal{R}}}_{\mathrm{SB}}({\mathcal{E}})$. The parameter n determines the rapidity of transition between the large-${\mathcal{R}}$ and small-${\mathcal{R}}$ regimes.

The same functional form might be adopted for g₁. Rather than make that choice, we first consider another possible constraint on g₁ and g₂.

The ${\mathcal{R}}$-directed flux coefficients that appear in the Fokker–Planck equation are given by Equations (6):

$$\begin{eqnarray}{D}_{{\mathcal{R}}} & = & -\langle {\rm{\Delta }}{\mathcal{R}}\rangle -\displaystyle \frac{5}{4{\mathcal{E}}}\langle {\rm{\Delta }}{\mathcal{E}}{\rm{\Delta }}{\mathcal{R}}\rangle \\ & & +\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {\mathcal{E}}}\langle {\rm{\Delta }}{\mathcal{E}}{\rm{\Delta }}{\mathcal{R}}\rangle +\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {\mathcal{R}}}\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle \end{eqnarray} \tag{ 26a }$$

$$\begin{eqnarray}&&\approx -\langle {\rm{\Delta }}{\mathcal{R}}\rangle +\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {\mathcal{R}}}\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle ,\end{eqnarray} \tag{ 26b }$$

$$\begin{eqnarray}&&{D}_{{\mathcal{R}}{\mathcal{R}}}=\displaystyle \frac{1}{2}\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle .\end{eqnarray} \tag{ 26c }$$

Since the ${\mathcal{R}}$-directed flux is

$$\begin{eqnarray}{\varphi }_{{\mathcal{R}}} & = & -{D}_{{\mathcal{R}}{\mathcal{E}}}\displaystyle \frac{\partial f}{\partial {\mathcal{E}}}-{D}_{{\mathcal{R}}{\mathcal{R}}}\displaystyle \frac{\partial f}{\partial {\mathcal{R}}}-{D}_{{\mathcal{R}}}f\\ & \approx & -{D}_{{\mathcal{R}}{\mathcal{R}}}\displaystyle \frac{\partial f}{\partial {\mathcal{R}}}-{D}_{{\mathcal{R}}}f,\end{eqnarray} \tag{ 27 }$$

it is reasonable to require that the diffusion coefficients in ${\mathcal{R}}$ satisfy

$$\begin{eqnarray}&&{D}_{{\mathcal{R}}}\to 0,\ \ \ {D}_{{\mathcal{R}}{\mathcal{R}}}\to 0\end{eqnarray} \tag{ 28 }$$

at the boundaries ${\mathcal{R}}=\{0,1\}$; in other words, that

$$\begin{eqnarray}&&\langle {\rm{\Delta }}{\mathcal{R}}\rangle =\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {\mathcal{R}}}\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle ,\ \langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle =0\end{eqnarray} \tag{ 29 }$$

at ${\mathcal{R}}=\{0,1\}$. Both the classical (Cohn–Kulsrud) and the resonant diffusion coefficients adopted here and in Papers I and II satisfy these conditions.

Suppose that a stronger condition is imposed: ${D}_{{\mathcal{R}}}\equiv 0$, $0\leqslant {\mathcal{R}}\leqslant 1$. In this case, the Fokker–Planck equation describing diffusion in angular momentum reduces to

$$\begin{eqnarray}\displaystyle \frac{\partial f}{\partial t} & = & -\displaystyle \frac{\partial {\varphi }_{{\mathcal{R}}}}{\partial {\mathcal{R}}}=\displaystyle \frac{\partial }{\partial {\mathcal{R}}}({D}_{{\mathcal{R}}{\mathcal{R}}}\displaystyle \frac{\partial f}{\partial {\mathcal{R}}})\\ & = & \displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {\mathcal{R}}}(\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle \displaystyle \frac{\partial f}{\partial {\mathcal{R}}})\end{eqnarray} \tag{ 30 }$$

which has a steady-state (zero-flux) solution $f({\mathcal{R}})=\mathrm{constant}$ regardless of the functional form of $\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle$. It could be argued that $f({\mathcal{R}})\;=$ constant is a reasonable form for a time-independent f, since it corresponds to an isotropic, or "maximum entropy," state. The resonant diffusion coefficients adopted in Papers I and II satisfy this stronger condition. Setting ${D}_{{\mathcal{R}}}\equiv 0$ also implies zero drift, i.e., zero flux in ${\mathcal{R}}$ in the absence of a gradient. For this reason, ${D}_{{\mathcal{R}}}\equiv 0$ will henceforth be called a "zero-drift" condition.

Returning now to the anomalous diffusion coefficients, we ask: what functional form for g₁ is required for zero drift? That is:

$$\begin{eqnarray}&&0={D}_{{\mathcal{R}}}=-{g}_{1}({\mathcal{R}})\langle {\rm{\Delta }}{\mathcal{R}}{\rangle }_{\mathrm{RR}}+\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {\mathcal{R}}}[{g}_{2}({\mathcal{R}})\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{RR}}].\end{eqnarray} \tag{ 31 }$$

Assuming that g₂ is given by Equation (25), the result is

$$\begin{eqnarray}{g}_{1}({\mathcal{R}}) & = & {\{1+{[\displaystyle \frac{{{\mathcal{R}}}_{2}({\mathcal{E}})}{{\mathcal{R}}}]}^{n}\}}^{-2/n}\\ & & +\displaystyle \frac{2(1-{\mathcal{R}})}{1-2{\mathcal{R}}}{(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{\mathcal{R}}})}^{n}{\{1+{[\displaystyle \frac{{{\mathcal{R}}}_{2}({\mathcal{E}})}{{\mathcal{R}}}]}^{n}\}}^{-2/n-1}.\end{eqnarray} \tag{ 32 }$$

Since the zero-drift argument is one of plausibility only, we will consider a slightly more general expression for g₁ that includes "zero drift" as a special case:

$$\begin{eqnarray}{g}_{1}({\mathcal{R}}) & = & {\{1+{[\displaystyle \frac{{{\mathcal{R}}}_{1}({\mathcal{E}})}{{\mathcal{R}}}]}^{n}\}}^{-2/n}\\ & & +\displaystyle \frac{2(1-{\mathcal{R}})}{1-2{\mathcal{R}}}{(\displaystyle \frac{{{\mathcal{R}}}_{1}}{{\mathcal{R}}})}^{n}{\{1+{[\displaystyle \frac{{{\mathcal{R}}}_{1}({\mathcal{E}})}{{\mathcal{R}}}]}^{n}\}}^{-2/n-1}.\end{eqnarray} \tag{ 33 }$$

The only difference between Equations (32) and (33) is the introduction of a second parameter, ${{\mathcal{R}}}_{1}$, in place of ${{\mathcal{R}}}_{2}$. This generalization still satisfies the zero-flux condition (28) at ${\mathcal{R}}=0$, regardless of ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$. At ${\mathcal{R}}=1$, g₁ and g₂ are both very close to one (especially since n will be chosen to be large) so that the resonant diffusion coefficients are recovered and the zero-flux condition is satisfied, again for any choice of ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$.

Correspondence of these expressions with the diffusion coefficients of Equation (22) at small ${\mathcal{R}}$ would require

$$\begin{eqnarray}{{\mathcal{R}}}_{1}^{2}({\mathcal{E}}) & = & (6/5)\tau A({\mathcal{E}}),\\ {{\mathcal{R}}}_{2}^{2}({\mathcal{E}}) & = & \tau A({\mathcal{E}}),\end{eqnarray} \tag{ 34 }$$

where

$$\begin{eqnarray}&&\tau A({\mathcal{E}})=4{\alpha }_{s}^{2}{(\displaystyle \frac{{t}_{\mathrm{coh}}}{P})}^{2}{(\displaystyle \frac{{r}_{g}}{a})}^{2}.\end{eqnarray} \tag{ 35 }$$

To within factors of order unity, these relations imply ${{\mathcal{R}}}_{1}\approx {{\mathcal{R}}}_{2}\approx {{\mathcal{R}}}_{\mathrm{SB}}^{({ii})}$ (Equation (18)).

It is shown in the Appendix that at small ${\mathcal{R}}$, these choices for g₁ and g₂ imply

$$\begin{eqnarray}&&{D}_{{\mathcal{R}}}\to 6\lambda A({\mathcal{E}}){(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{2},\ \ \ \lambda \equiv 1-\displaystyle \frac{{{\mathcal{R}}}_{2}^{2}}{{{\mathcal{R}}}_{1}^{2}}\end{eqnarray} \tag{ 36 }$$

so that the direction of the drift is determined by the relative sizes of ${{\mathcal{R}}}_{1}$ and ${{\mathcal{R}}}_{2}$, as follows:

$$\begin{eqnarray}&&{{\mathcal{R}}}_{1}\lt {{\mathcal{R}}}_{2}\ \to \ {D}_{{\mathcal{R}}}\lt 0\ \to \ {\varphi }_{{\mathcal{R}}}\gt 0\end{eqnarray} \tag{ 37a }$$

$$\begin{eqnarray}&&{{\mathcal{R}}}_{1}\gt {{\mathcal{R}}}_{2}\ \to \ {D}_{{\mathcal{R}}}\gt 0\ \to \ {\varphi }_{{\mathcal{R}}}\lt 0\end{eqnarray} \tag{ 37b }$$

where the expressions for ${\varphi }_{{\mathcal{R}}}$ assume $f({\mathcal{R}})=\mathrm{constant}$. Furthermore the form of both the steady-state $f({\mathcal{R}})$ and the steady-state flux (assuming the presence of a sink, i.e., that $f({\mathcal{R}})=0$ at ${\mathcal{R}}\leqslant {{\mathcal{R}}}_{0}$) depend sensitively on λ, as shown in Figures 10 and 11 from the Appendix. In the zero-drift ($\lambda =0$) case, the steady-state flux is reduced by a factor

$$\begin{eqnarray}&&\eta \approx 2{(\displaystyle \frac{{{\mathcal{R}}}_{0}}{{{\mathcal{R}}}_{2}})}^{2}\mathrm{log}(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{{\mathcal{R}}}_{0}})\end{eqnarray} \tag{ 38 }$$

compared with the flux that would be obtained in the absence of anomalous relaxation (note that an empty loss cone has been assumed). For nonzero λ, the reduction factor can be greater or smaller than this, tending toward a maximum value of one for large ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$ (Figure 11).

Bar-Or & Alexander (2014) suggested a different functional form for $\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle$ in the anomalous regime:

$$\begin{eqnarray}&&\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle \propto \mathrm{exp}(-\displaystyle \frac{1}{\pi }\displaystyle \frac{{{\mathcal{R}}}_{\mathrm{SB}}^{2}}{{{\mathcal{R}}}^{2}}),\end{eqnarray} \tag{ 39 }$$

an exponential cut-off toward small ${\mathcal{R}}$. While there does not seem to be strong support for this functional form in any published numerical simulations, we consider it here for completeness, and because it serves to illustrate how sensitively the evolution of f depends on the form assumed for the anomalous diffusion coefficients.

Proceeding as above, we write

$$\begin{eqnarray}\langle {\rm{\Delta }}{\mathcal{R}}\rangle & = & \langle {\rm{\Delta }}{\mathcal{R}}{\rangle }_{\mathrm{CK}}+{h}_{1}({\mathcal{E}},{\mathcal{R}})\langle {\rm{\Delta }}{\mathcal{R}}{\rangle }_{\mathrm{RR}},\\ \langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle & = & \langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{CK}}+{h}_{2}({\mathcal{E}},{\mathcal{R}})\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{RR}}.\end{eqnarray} \tag{ 40 }$$

Suppose

$$\begin{eqnarray}&&{h}_{2}({\mathcal{E}},{\mathcal{R}})=\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}})\end{eqnarray} \tag{ 41 }$$

where ${{\mathcal{R}}}_{4}({\mathcal{E}})\approx {{\mathcal{R}}}_{\mathrm{SB}}({\mathcal{E}})$. Since ${{\mathcal{R}}}_{4}\ll 1$, ${h}_{2}({\mathcal{E}},1)\approx 1$. The zero-drift condition would imply

$$\begin{eqnarray*}&&\langle {\rm{\Delta }}{\mathcal{R}}\rangle =2A({\mathcal{E}})[1-2{\mathcal{R}}+2(1-{\mathcal{R}}){(\displaystyle \frac{{{\mathcal{R}}}_{4}}{{\mathcal{R}}})}^{2}]\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}}).\end{eqnarray*}$$

We again generalize this expression by defining a second parameter, ${{\mathcal{R}}}_{3}$, and writing

$$\begin{eqnarray}&&{h}_{1}({\mathcal{R}})=[1+\displaystyle \frac{2(1-{\mathcal{R}})}{1-2{\mathcal{R}}}{(\displaystyle \frac{{{\mathcal{R}}}_{3}}{{\mathcal{R}}})}^{2}]\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{3}^{2}}{{{\mathcal{R}}}^{2}}).\end{eqnarray} \tag{ 42 }$$

At small ${\mathcal{R}}$, these expressions imply a flux coefficient

$$\begin{eqnarray}{D}_{{\mathcal{R}}} & \to & 4A({\mathcal{E}}){(\displaystyle \frac{{{\mathcal{R}}}_{4}}{{\mathcal{R}}})}^{2}\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}})\\ & & \times [1-\displaystyle \frac{{{\mathcal{R}}}_{3}^{2}}{{{\mathcal{R}}}_{4}^{2}}\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{3}^{2}-{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}})]\end{eqnarray} \tag{ 43 }$$

showing that in this case, as in the power-law case, the sign of ${D}_{{\mathcal{R}}}$ is determined by the relative sizes of ${{\mathcal{R}}}_{3}$ and ${{\mathcal{R}}}_{4}$. Figures 10 and 11 illustrate the properties of the steady-state solutions $f({\mathcal{R}})$. The behavior of f at small ${\mathcal{R}}$ depends very sensitively on ${{\mathcal{R}}}_{3}/{{\mathcal{R}}}_{4}$. The reduction in the steady-state flux, for ${{\mathcal{R}}}_{3}={{\mathcal{R}}}_{4}$, is

$$\begin{eqnarray}&&\eta \approx 2\displaystyle \frac{{{\mathcal{R}}}_{4}}{{{\mathcal{R}}}_{0}}\mathrm{log}(\displaystyle \frac{{{\mathcal{R}}}_{4}}{{{\mathcal{R}}}_{0}}){e}^{-{{\mathcal{R}}}_{4}/{{\mathcal{R}}}_{0}}.\end{eqnarray} \tag{ 44 }$$

Based on the analysis presented above and in the Appendix, many properties of the steady-state solutions are expected to depend sensitively on the functional forms of the diffusion coefficients in the anomalous regime, and in particular on the ratios ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$ or ${{\mathcal{R}}}_{3}/{{\mathcal{R}}}_{4}$.

In Paper I, the functional forms of the diffusion coefficients in the resonant regime were successfully constrained by comparison with the numerical results of Hamers et al. (2014). Those authors used a test-particle integrator to extract values of the angular momentum diffusion coefficients for stars orbiting near a SBH in nuclei with $n(r)\propto {r}^{-2}$ and ${r}^{-1}$.

The Hamers et al. integrations included post-Newtonian terms in the equations of motion, both for the field and test stars, and the suppression of angular momentum diffusion below the SB was clearly seen. Figure 1 provides an illustration: it shows the first- and second-order diffusion coefficients for stars in a single energy bin, in integrations of a model with $n(r)\propto {r}^{-2}$.¹ Overplotted are analytic diffusion coefficients from the two families considered above: power-law (Equations (66a), (67)) and exponential (Equations (66a), (88)).

Zoom In Zoom Out Reset image size

These figures, and similar ones made for stars at other energies, motivate the following conclusions (some of which were presented already by Hamers et al.).

• 1.  
  An exponential dependence of the diffusion coefficients on ${\mathcal{R}}$ in the anomalous regime (Bar-Or & Alexander 2014) is ruled out.
• 2.  
  The power-law dependence of $\langle {\rm{\Delta }}{\mathcal{R}}\rangle$ and $\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle$ on ${\mathcal{R}}$ predicted in Paper I, and reproduced here in Equations (22), is confirmed, particularly in the case of the second-order coefficient for which the noise is smallest.
• 3.  
  The value of ${\mathcal{R}}$ that defines the transition between the resonant and anomalous regimes, called here ${{\mathcal{R}}}_{2}$, is well predicted by Equation (18).

An attempt was made to find the best-fit value(s) of ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$ by searching for parameter values that optimized the fits to data like those in Figure 1. Unfortunately, the results so obtained were found not to be robust. This was due in part to the greater noise associated with data below the SB, particularly in the case of the first-order coefficient. In addition, the much greater variation in the amplitudes of $\langle {\rm{\Delta }}{\mathcal{R}}\rangle$ and $\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle$ at a single energy meant that the best-fit solution depended sensitively on the relative weighting of the data at different values of ${\mathcal{R}}$. In the end, all that could be concluded was that the first-order diffusion coefficients are consistent with ${{\mathcal{R}}}_{1}={{\mathcal{R}}}_{2}$, the "zero-drift" condition, but that values of ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$ moderately greater or less than one could not be excluded using these data.

The fact that the steady-state form of $f({\mathcal{R}})$, and the loss-cone flux, depends sensitively on ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$ suggests a second way to constrain the anomalous diffusion coefficients: insert them into the Fokker–Planck equation and integrate forward from initial conditions like the ones used by Merritt et al. (2011) in their small-N simulations. As discussed in Hamers et al. the value of N in those simulations was too small to allow direct extraction of the diffusion coefficients. However, the time-averaged, or integrated, properties of the N-body models were reasonably well determined, particularly given that multiple (∼8) realizations of the same initial conditions were integrated, allowing the variance in the results to be reduced by averaging.

The Merritt et al. (2011) initial conditions consisted of 50 stars, of mass ${m}_{\star }=50{M}_{\odot }$, distributed as $n(r)\propto {r}^{-2}$ around a SBH of mass ${M}_{\bullet }={10}^{6}{M}_{\odot }$. The initial distribution was truncated for orbits with semimajor axes above 10 mpc and below 0.1 mpc. Integrations were carried out both with, and without, post-Newtonian terms in the equations of motion, up to order 2.5 PN. Capture of stars by the SBH was allowed to occur whenever the orbital periapsis fell below $8{r}_{g}=8{{GM}}_{\bullet }/{c}^{2}$. Each realization of the initial conditions was integrated for a time corresponding to $2\times {10}^{6}$ year (with PN terms) and 10⁷ year (without PN terms). The average capture rate in the relativistic integrations was about one event per 10⁶ year; in the absence of the PN terms, mean capture rates were about a factor 20 higher.

There is no ambiguity in representing the Merritt et al. N-body initial conditions as a smooth $f({\mathcal{E}},{\mathcal{R}})$. However, the Fokker–Planck algorithm has a number of parameters that must be specified, in addition to those that define the anomalous diffusion coefficients (n and ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$, or ${{\mathcal{R}}}_{3}/{{\mathcal{R}}}_{4}$). Those parameters include

$$\begin{eqnarray}&&{{\mathcal{E}}}_{\mathrm{min}},\ \ \ \mathrm{ln}{\rm{\Lambda }},\ \ \ \displaystyle \frac{{r}_{\mathrm{lc}}}{{r}_{g}}\ \ \ \ {\rm{\Delta }}t,\ \ \ {N}_{X},\ \ \ {N}_{Z}.\end{eqnarray} \tag{ 45 }$$

${{\mathcal{E}}}_{\mathrm{min}}$ is the binding energy at the edge of the (${\mathcal{E}},{\mathcal{R}}$) grid; it should be small enough that few stars diffuse to ${\mathcal{E}}\lt {{\mathcal{E}}}_{\mathrm{min}}$ during the course of the integration. The value ${10}^{-8}{c}^{2}$ was chosen, which is the energy of an orbit with semimajor axis $a=5\times {10}^{7}{r}_{g}\;\approx$ 2.5 pc, or ∼250 times the maximum a-value of the initial conditions. The number of grid points was ${N}_{X}={N}_{Z}=64$. The quantity $\mathrm{ln}{\rm{\Lambda }}$ was set to 15 in most of the integrations, except for one set in which smaller and larger values (from 11 to 19) were tried. The Coulomb logarithm only appears in the expressions for the classical diffusion coefficients; since evolution of these models is dominated by resonant relaxation, the results are expected to be weakly dependent on $\mathrm{ln}{\rm{\Lambda }}$ and this was found to be the case. The integration time step, ${\rm{\Delta }}t$, was set to 2000 years, i.e., 10³ steps per integration.

A natural choice for the parameter ${r}_{\mathrm{lc}}/{r}_{g}$ in the Fokker–Planck integrations would be 8, the same value assumed in the N-body integrations. In the case of "plunges"—captures that occur without significant energy loss due to gravitational radiation—this would be the correct choice. However, some of the N-body capture events were "EMRIs," for which capture was preceded by significant energy loss due to the 2.5PN terms. No such loss terms are included in the Fokker–Planck integrations described here. Roughly speaking, the effect of the 2.5PN terms is to shift the location of the loss cone toward larger ${\mathcal{R}}$ (i.e., larger orbital periapsis) at each ${\mathcal{E}}$ (see e.g., Figure 5 of Merritt et al. 2011). To evaluate the effect on the Fokker–Planck results of ignoring the 2.5PN terms, a set of integrations was carried out setting ${r}_{\mathrm{lc}}/{r}_{g}=32$, four times its value in the N-body integrations.

Figure 2 compares time-averaged loss rates in the N-body and Fokker–Planck models, defined as the total number of stars lost until time t divided by t. The left panel shows results using the power-law forms of the anomalous diffusion coefficients, Equations (66a), (67), with n = 8; the right panel shows results using the exponential forms, Equations (66a), (88). Each panel shows results for several values of ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$ (power-law) or ${{\mathcal{R}}}_{3}/{{\mathcal{R}}}_{4}$ (exponential), as specified in the caption. For each value of this ratio, three integrations were carried out, varying the way in which ${{\mathcal{R}}}_{2}$ or ${{\mathcal{R}}}_{4}$ were related to ${{\mathcal{R}}}_{\mathrm{SB}}$. One of the three integrations (shown by the solid curves) equated ${{\mathcal{R}}}_{2}$ or ${{\mathcal{R}}}_{4}$ with ${{\mathcal{R}}}_{\mathrm{SB}}^{({ii})}$, Equation (18). The other two integrations adopted larger or smaller values: by a factor two or one-half (in the power-law models), or by a factor $\sqrt{\pi }$ or $1/\sqrt{\pi }$ (exponential). These additional integrations are shown by the dashed–dotted curves in Figure 2. Larger (smaller) values of ${{\mathcal{R}}}_{2}$ or ${{\mathcal{R}}}_{4}$ generally resulted in smaller (larger) loss rates, at least at early times.

Zoom In Zoom Out Reset image size

Figure 2 suggests that both functional forms of the anomalous diffusion coefficients are able to reproduce the N-body capture rates, as long as ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$ or ${{\mathcal{R}}}_{3}/{{\mathcal{R}}}_{4}$ are not too different from one. The best correspondence is achieved, in both cases, when this ratio is slightly less than one, and this result remains unchanged even when the values of ${{\mathcal{R}}}_{2}({\mathcal{E}})$ or ${{\mathcal{R}}}_{4}({\mathcal{E}})$ are substantially modified. Recall that ${{\mathcal{R}}}_{1}\lt {{\mathcal{R}}}_{2}$ or ${{\mathcal{R}}}_{3}\lt {{\mathcal{R}}}_{4}$ imply ${D}_{{\mathcal{R}}}\lt 0$, i.e., ${\varphi }_{{\mathcal{R}}}\gt 0$, i.e., a drift toward larger ${\mathcal{R}}$ (Equation (37)).

Changing the parameter n in the power-law diffusion coefficients from n = 8 to n = 32 had almost no discernible effect on the loss rates. Varying $\mathrm{ln}{\rm{\Lambda }}$ or ${r}_{\mathrm{lc}}/{r}_{g}$ in the amounts described above did result in noticeable changes, but by amounts comparable with the ranges shown in Figure 2 due to variations in the definition of ${{\mathcal{R}}}_{2}$ or ${{\mathcal{R}}}_{4}$. In every set of integrations, correspondence with the N-body loss rates was best for ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}\lesssim 1$ or ${{\mathcal{R}}}_{3}/{{\mathcal{R}}}_{4}\lesssim 1$.

Figure 3 makes another comparison between Fokker–Planck and N-body models. Plotted there are time-averaged angular momentum distributions at a single energy. These are displayed as ${dN}/{dX}$, where X is defined as in Figure 11 of Merritt et al. (2011):

$$\begin{eqnarray}&&X\equiv {\mathrm{log}}_{10}(1-e)\end{eqnarray} \tag{ 46 }$$

with $e=\sqrt{1-{\mathcal{R}}}$ the orbital eccentricity. Figure 3 implies that values of ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$ or ${{\mathcal{R}}}_{3}/{{\mathcal{R}}}_{4}$ of unity ("zero drift") or greater can be securely ruled out—consistent with Figure 2. In the case of the power-law forms of the diffusion coefficients (upper panels), the best correspondence with the N-body results seems to occur for

$$\begin{eqnarray}&&{{\mathcal{R}}}_{2}\approx 2{{\mathcal{R}}}_{\mathrm{SB}}^{({ii})},\ \ \ 0.8\lesssim {{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}\lesssim 0.9.\end{eqnarray} \tag{ 47 }$$

In the case of the exponential forms of the diffusion coefficients (lower panels), correspondence with the N-body results seems to require

$$\begin{eqnarray}&&{{\mathcal{R}}}_{2}\approx \sqrt{\pi }\ {{\mathcal{R}}}_{\mathrm{SB}}^{({ii})},\ \ \ 0.9\lesssim {{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}\lesssim 0.95.\end{eqnarray} \tag{ 48 }$$

Once again, the best correspondence is achieved with diffusion coefficients that imply a non-zero drift, in the direction of increasing ${\mathcal{R}}$, in the anomalous regime.

Zoom In Zoom Out Reset image size

Figure 4 shows representative plots of $f({\mathcal{R}};{\mathcal{E}})$ at one ${\mathcal{E}}$ from the Fokker–Planck models at the final time step (roughly $2.5\times {10}^{6}$ year). Dashed curves show models having parameters similar to those found to correspond best to the N-body results. These solutions always exhibit a rapid drop in the steady-state $f({\mathcal{R}})$ below the SB. That drop would be even steeper in the absence of classical relaxation, which dominates the evolution in ${\mathcal{R}}$ at small ${\mathcal{R}}$ (the region left of the vertical dotted lines in the figure), thus maintaining a relatively high diffusion rate at small ${\mathcal{R}}$.

Zoom In Zoom Out Reset image size

One final argument can be made in support of diffusion coefficients that satisfy ${D}_{{\mathcal{R}}}\lt 0$. As described in Antonini & Merritt (2013), stars orbiting near the SB are observed to exhibit a "buoyancy" phenomenon: should they cross the barrier from below (${\mathcal{R}}\lt {{\mathcal{R}}}_{\mathrm{SB}}$) to above (${\mathcal{R}}\gt {{\mathcal{R}}}_{\mathrm{SB}}$), they tend to remain above, and vice- versa. This behavior is consistent with a drift toward larger ${\mathcal{R}}$, as implied by ${D}_{{\mathcal{R}}}\lt 0$.

In the regime of anomalous relaxation (${\mathcal{R}}\lt {{\mathcal{R}}}_{\mathrm{SB}}({\mathcal{E}})$), timescales for angular momentum diffusion become long for small ${\mathcal{R}}$. One consequence is that classical relaxation, which (by assumption) is not affected by precession, can once again dominate the evolution in angular momentum (Hamers et al. 2014).

We estimate the value of ${\mathcal{R}}$ at which this occurs by equating the first and second terms on the right hand sides of Equation (24).

We simplify the expressions by using the limiting forms of the diffusion coefficients as ${\mathcal{R}}\to 0$.

For the classical diffusion coefficients, Equations (11) and (12) from Hamers et al. (2014), together with the transformation Equation (32) from Paper I, yield

$$\begin{eqnarray}&&\langle {\rm{\Delta }}{\mathcal{R}}{\rangle }_{\mathrm{CK}}\to {A}^{H}({\mathcal{E}}),\ \ \ \langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{CK}}\to 2{\mathcal{R}}{A}^{H}({\mathcal{E}}),\end{eqnarray} \tag{ 49a }$$

$$\begin{eqnarray}&&{A}^{H}(a)=\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{{C}_{\mathrm{NRR}}(\gamma )}{(\displaystyle \frac{{m}_{\star }}{{M}_{\bullet }})}^{2}\displaystyle \frac{N(a)}{P(a)}\end{eqnarray} \tag{ 49b }$$

where the symbol "$\to$" denotes the limit of small ${\mathcal{R}}$ and the function ${C}_{\mathrm{NRR}}(\gamma )$ is given in Appendix B of Hamers et al. (2014); their calculation assumes $\rho (r)\propto {r}^{-\gamma }$.

For the anomalous diffusion coefficients, Equations (24), (25) and (15) give the limiting forms

$$\begin{eqnarray}\langle {\rm{\Delta }}{\mathcal{R}}{\rangle }_{\mathrm{AR}} & \to & 6A({\mathcal{E}}){(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{1}})}^{2},\\ \langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{AR}} & \to & 4A({\mathcal{E}}){\mathcal{R}}{(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{2}\end{eqnarray} \tag{ 50 }$$

in the power-law case, with $A(a)$ defined in Equation (15). Equating $\langle {\rm{\Delta }}{\mathcal{R}}{\rangle }_{\mathrm{CK}}$ with $\langle {\rm{\Delta }}{\mathcal{R}}{\rangle }_{\mathrm{AR}}$ or $\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{CK}}$ with $\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{AR}}$ yields

$$\begin{eqnarray}&&\displaystyle \frac{{{\mathcal{R}}}^{2}}{{{\mathcal{R}}}_{1}^{2}}=\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{6\;{C}_{\mathrm{NR}}{\alpha }^{2}}\displaystyle \frac{P}{{t}_{\mathrm{coh}}}\ \ \ \ (1\mathrm{st}\ \mathrm{order}),\end{eqnarray} \tag{ 51a }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{\mathcal{R}}}^{2}}{{{\mathcal{R}}}_{2}^{2}}=\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{2\;{C}_{\mathrm{NR}}{\alpha }^{2}}\displaystyle \frac{P}{{t}_{\mathrm{coh}}}\ \ \ \ (2\mathrm{nd}\ \mathrm{order}).\end{eqnarray} \tag{ 51b }$$

Replacing ${{\mathcal{R}}}_{1}$ and ${{\mathcal{R}}}_{2}$ by ${{\mathcal{R}}}_{\mathrm{SB}}^{({ii})}$ in Equation (51a) and setting $\alpha =1.6$ yields

$$\begin{eqnarray}&&{{\mathcal{R}}}^{2}=(1.0,3.1)\ \displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{{C}_{\mathrm{NRR}}}{(\displaystyle \frac{{r}_{g}}{a})}^{2}\displaystyle \frac{{t}_{\mathrm{coh}}}{P}\end{eqnarray} \tag{ 52 }$$

where the constants in parentheses refer to first- and second-order diffusion coefficients respectively. Equation (52) is similar to Equation (23a) of Hamers et al. (2014).

Equations (51a) (51b) were used to plot the vertical dotted lines in the left-hand panels of Figures 1 and 4 (note that these two figures refer to different mass models). In Figure 1, the line predicts reasonably well the value of ${\mathcal{R}}$ at which the data begin to deviate from the fitted curves. In the case of Figure 4, the effects of classical relaxation can be seen in the steady-state $f({\mathcal{R}})$, which drops more gradually to zero below the SB than it would if only anomalous relaxation were acting (Figure 10).

If the exponential forms of the anomalous diffusion coefficients are correct, then the value of ${\mathcal{R}}$ at which $\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{CK}}=\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}{\rangle }_{\mathrm{AR}}$ becomes

$$\begin{eqnarray}&&\displaystyle \frac{{{\mathcal{R}}}^{2}}{{{\mathcal{R}}}_{4}^{2}}={[\mathrm{log}(\displaystyle \frac{2{\alpha }^{2}{C}_{\mathrm{NR}}}{\mathrm{ln}{\rm{\Lambda }}}\displaystyle \frac{{t}_{\mathrm{coh}}}{P})]}^{-1}.\end{eqnarray} \tag{ 53 }$$

This value for ${\mathcal{R}}$ is plotted, in addition to the value given by Equation (52), on the lower right-hand panel of Figure 1. Note that under this hypothesis, the range in ${\mathcal{R}}$ over which anomalous relaxation would be relevant would become very small; in effect, classical relaxation would dominate the evolution everywhere below (and sometimes even above) the SB. We reiterate that there is no support for the exponential form of the anomalous diffusion coefficients in any published numerical simulations.

In Paper II, the formation of cores due to resonant relaxation was discussed. Equation (38) of that paper gave an estimate of the radius of the core so formed:

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{\mathrm{core}}}{{r}_{{\rm{m}}}}\approx 0.028{(\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{15})}^{-4/5}.\end{eqnarray} \tag{ 58 }$$

Here, r_m is the gravitational influence radius of the SBH, defined as the radius of a sphere containing a mass in stars of $2{M}_{\bullet }$.

As described in this paper, at least in the case of the Milky Way, the inclusion of anomalous relaxation tends to counteract core formation by causing stars to accumulate near the SB.

Here we consider the generality of that result. Suppose that nuclei of galaxies fainter than the Milky Way have SBHs that satisfy the ${M}_{\bullet }-\sigma$ relation:

$$\begin{eqnarray*}&&\displaystyle \frac{{M}_{\bullet }}{{10}^{6}{M}_{\odot }}\approx 5.76{(\displaystyle \frac{\sigma }{100\;\mathrm{km}\ {{\rm{s}}}^{-1}})}^{4.86}\end{eqnarray*}$$

(Merritt 2013, Equation (2.33)) and that their nuclear densities are close to the Bahcall–Wolf form, $\rho (r)\propto {r}^{-7/4}$. The latter is not too different from $\rho \propto {r}^{-2}$, for which ${r}_{{\rm{m}}}\approx {r}_{h}\equiv {{GM}}_{\bullet }/{\sigma }^{2}$. Under these assumptions,

$$\begin{eqnarray}&&{r}_{{\rm{m}}}\approx 0.88{(\displaystyle \frac{{M}_{\bullet }}{{10}^{6}{M}_{\odot }})}^{0.59}\mathrm{pc},\ \ \ {M}_{\bullet }\lesssim 4\times {10}^{6}{M}_{\odot }.\end{eqnarray} \tag{ 59 }$$

We first verify the assumptions that led to Equation (58). That equation was derived assuming that ${t}_{\mathrm{coh}}={t}_{\mathrm{coh},{\rm{M}}}$ (Equation (16)). The radius at which ${t}_{\mathrm{coh},{\rm{M}}}={t}_{\mathrm{coh},{\rm{S}}}$, for a nucleus with $\rho (r)\propto {r}^{-7/4}$, is given by Equation (30) from Paper II:

$$\begin{eqnarray}\displaystyle \frac{{a}_{\mathrm{coh}}}{{r}_{{\rm{m}}}} & \approx & 1.2\times {10}^{-3}{(\displaystyle \frac{{M}_{\bullet }}{{10}^{6}{M}_{\odot }})}^{4/9}{(\displaystyle \frac{{r}_{{\rm{m}}}}{1\mathrm{pc}})}^{-4/9}\\ & \approx & 1.3\times {10}^{-3}{(\displaystyle \frac{{M}_{\bullet }}{{10}^{6}{M}_{\odot }})}^{0.18}\end{eqnarray} \tag{ 60 }$$

suggesting that indeed ${t}_{\mathrm{coh}}\sim {t}_{\mathrm{coh},{\rm{M}}}$ at $a={a}_{\mathrm{core}}$ in the galaxies of interest.

Next we ask how a_core compares with the radii that define the SB. Equation (20a) gave approximate limits on the extent of the SB in a $\rho \propto {r}^{-7/4}$ nucleus, which we recast here as

$$\begin{eqnarray}\displaystyle \frac{{a}_{\mathrm{min}}}{{r}_{{\rm{m}}}} & \approx & {(\displaystyle \frac{{M}_{\bullet }}{2{m}_{\star }})}^{4/13}{(\displaystyle \frac{{r}_{g}}{{r}_{{\rm{m}}}})}^{8/13},\\ \displaystyle \frac{{a}_{\mathrm{max}}}{{r}_{{\rm{m}}}} & \approx & {(4{\rm{\Theta }})}^{-4/9}{(\displaystyle \frac{{M}_{\bullet }}{{m}_{\star }})}^{4/9}{(\displaystyle \frac{{r}_{g}}{{r}_{{\rm{m}}}})}^{4/9}.\end{eqnarray} \tag{ 61 }$$

Comparing a_min and a_max to a_core:

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{\mathrm{min}}}{{a}_{\mathrm{core}}}\approx 29{(\displaystyle \frac{{M}_{\bullet }}{{m}_{\star }})}^{4/13}{(\displaystyle \frac{{r}_{g}}{{r}_{{\rm{m}}}})}^{8/13}{(\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{15})}^{4/5},\end{eqnarray} \tag{ 62a }$$

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{\mathrm{max}}}{{a}_{\mathrm{core}}}\approx 5.8{(\displaystyle \frac{{\rm{\Theta }}}{15})}^{-4/9}{(\displaystyle \frac{{M}_{\bullet }}{{m}_{\star }})}^{4/9}{(\displaystyle \frac{{r}_{g}}{{r}_{{\rm{m}}}})}^{4/9}{(\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{15})}^{4/5}.\end{eqnarray} \tag{ 62b }$$

Using ${r}_{g}\approx 4.78\times {10}^{-8}({M}_{\bullet }/{10}^{6}{M}_{\odot })\;\mathrm{pc}$, these become

$$\begin{eqnarray}\displaystyle \frac{{a}_{\mathrm{min}}}{{a}_{\mathrm{core}}} & \approx & 0.064{(\displaystyle \frac{{M}_{\bullet }}{{10}^{6}{M}_{\odot }})}^{12/13}{(\displaystyle \frac{{m}_{\star }}{{M}_{\odot }})}^{-4/13}\\ & & \times {(\displaystyle \frac{{r}_{{\rm{m}}}}{1\ \mathrm{pc}})}^{-8/13}{(\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{15})}^{4/5},\end{eqnarray} \tag{ 63a }$$

$$\begin{eqnarray}\displaystyle \frac{{a}_{\mathrm{max}}}{{a}_{\mathrm{core}}} & \approx & 1.5{(\displaystyle \frac{{\rm{\Theta }}}{15})}^{-4/9}{(\displaystyle \frac{{M}_{\bullet }}{{10}^{6}{M}_{\odot }})}^{8/9}\\ & & \times {(\displaystyle \frac{{m}_{\star }}{{M}_{\odot }})}^{-4/9}{(\displaystyle \frac{{r}_{{\rm{m}}}}{1\ \mathrm{pc}})}^{-4/9}{(\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{15})}^{4/5}.\end{eqnarray} \tag{ 63b }$$

Setting ${a}_{\mathrm{max}}\lt {a}_{\mathrm{core}}$ implies that anomalous relaxation is unlikely to affect the formation of the core due to resonant relaxation. This condition is:

$$\begin{eqnarray}&&{r}_{{\rm{m}}}\gtrsim 2.5{(\displaystyle \frac{{\rm{\Theta }}}{15})}^{-1}{(\displaystyle \frac{{M}_{\bullet }}{{10}^{6}{M}_{\odot }})}^{2}{(\displaystyle \frac{{m}_{\star }}{{M}_{\odot }})}^{-1}{(\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{15})}^{9/5}\mathrm{pc}.\end{eqnarray} \tag{ 64 }$$

In the case of the Milky Way (${M}_{\bullet }\approx 4\times {10}^{6}{M}_{\odot }$), satisfying this condition for ${m}_{\star }={M}_{\odot }$ would require ${r}_{{\rm{m}}}\gtrsim 40\ \mathrm{pc}$—about ten times larger than the value inferred from stellar kinematics. This is consistent with Figure 8, which showed that anomalous relaxation inhibits the formation of a core for all reasonable values of r_m. In the case of galaxies with central black holes less massive than the Milky Way's, Equations (64) and (59) allow us to write the condition ${a}_{\mathrm{max}}\lt {a}_{\mathrm{core}}$ as

$$\begin{eqnarray}&&{M}_{\bullet }\lesssim 5\times {10}^{5}{(\displaystyle \frac{{\rm{\Theta }}}{15})}^{0.48}{(\displaystyle \frac{{m}_{\star }}{{M}_{\odot }})}^{0.48}{(\displaystyle \frac{\mathrm{ln}{\rm{\Lambda }}}{15})}^{-0.86}{M}_{\odot }.\end{eqnarray} \tag{ 65 }$$

Thus, the core formed by resonant relaxation is expected to become progressively more prominent as ${M}_{\bullet }$ is reduced below its value in the Milky Way. This fact is likely to be important in determining the rate of formation of gravitational-wave sources, particularly since theoretical estimates often focus on galaxies with ${M}_{\bullet }\lesssim {10}^{6}{M}_{\odot }$. Estimating this rate will be the topic of upcoming papers in this series.

Until recently, discussions of gravitational encounters near a SBH have usually been presented in terms of diffusion timescales, that is, in terms of second-order diffusion coefficients like $\langle {({\rm{\Delta }}L)}^{2}\rangle$ (Rauch & Tremaine 1996; Hopman & Alexander 2006; Gürkan and Hopman 2007; Eilon et al. 2009; Madigan et al. 2011). Hamers et al. (2014) were apparently the first to consider the forms of the first-order diffusion coefficients. Of course, both first- and second-order diffusion coefficients are essential when computing the evolution of f via the Fokker–Planck equation.

As shown here through a number of examples, the form of the steady-state $f({\mathcal{E}},{\mathcal{R}})$ near the SB can depend very sensitively on the relative amplitude of the first- and second-order coefficients in the anomalous-relaxation regime. A natural case to consider is that of "zero drift," in which the two diffusion coefficients imply a net flux in angular momentum that is zero when $f({\mathcal{R}})=\mathrm{constant}$ (Section 3.1). While it may be natural, this assumption is not compelled by any physical argument of which we are aware. Furthermore, as discussed in detail in Section 3.3, the highest-quality N-body simulations carried out to date of this regime (Merritt et al. 2011) seem to require anomalous diffusion coefficients that differ slightly, though significantly, from the "zero-drift" condition. A state of "positive drift" seems to better characterize the existing simulations.

Elucidation of the long-term effects of gravitational encounters in the Schwarzschild regime near a SBH will ultimately require a better specification of the anomalous diffusion coefficients. By far the best way to do this—at least in principle—is via direct N-body integrations, which impose the fewest approximations. In practice, integrations of the required accuracy become very time-consuming when $N\gtrsim {10}^{2}$. A major effort should be devoted to increasing the efficiency of the N-body integrators.

Identify w₁ and w₂ with g₁ and g₂ given respectively by Equations (33) and (25):

$$\begin{eqnarray*}{g}_{1}({\mathcal{E}},{\mathcal{R}}) & = & {\{1+{[\displaystyle \frac{{{\mathcal{R}}}_{1}({\mathcal{E}})}{{\mathcal{R}}}]}^{n}\}}^{-2/n}\\ & & +\displaystyle \frac{2(1-{\mathcal{R}})}{1-2{\mathcal{R}}}{(\displaystyle \frac{{{\mathcal{R}}}_{1}}{{\mathcal{R}}})}^{n}{\{1+{[\displaystyle \frac{{{\mathcal{R}}}_{1}({\mathcal{E}})}{{\mathcal{R}}}]}^{n}\}}^{-2/n-1},\end{eqnarray*} \tag{ 67 }$$

$$\begin{eqnarray*}&&{g}_{2}({\mathcal{E}},{\mathcal{R}})={\{1+{[\displaystyle \frac{{{\mathcal{R}}}_{2}({\mathcal{E}})}{{\mathcal{R}}}]}^{n}\}}^{-2/n}.\end{eqnarray*}$$

The flux coefficients, Equations (26), are

$$\begin{eqnarray}&&{D}_{{\mathcal{R}}}=-\langle {\rm{\Delta }}{\mathcal{R}}\rangle +\displaystyle \frac{1}{2}\displaystyle \frac{\partial }{\partial {\mathcal{R}}}\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle \end{eqnarray} \tag{ 68a }$$

$$\begin{eqnarray} & = & 2A({\mathcal{E}}){[1+{(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{\mathcal{R}}})}^{n}]}^{-(1+2/n)}[1-2{\mathcal{R}}+(3-4{\mathcal{R}}){(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{\mathcal{R}}})}^{n}]\\ & & -2A({\mathcal{E}}){[1+{(\displaystyle \frac{{{\mathcal{R}}}_{1}}{{\mathcal{R}}})}^{n}]}^{-(1+2/n)}\\ & & \times [1-2{\mathcal{R}}+(3-4{\mathcal{R}}){(\displaystyle \frac{{{\mathcal{R}}}_{1}}{{\mathcal{R}}})}^{n}]\end{eqnarray} \tag{ 68b }$$

$$\begin{eqnarray}&&{D}_{{\mathcal{R}}{\mathcal{R}}}=\displaystyle \frac{1}{2}\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle \end{eqnarray} \tag{ 68c }$$

$$\begin{eqnarray}&&=\;2A({\mathcal{E}}){\mathcal{R}}(1-{\mathcal{R}}){[1+{(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{\mathcal{R}}})}^{n}]}^{-2/n}.\end{eqnarray} \tag{ 68d }$$

It is easy to verify that in this case,

$$\begin{eqnarray}&&{D}_{{\mathcal{R}}}=-\langle {\rm{\Delta }}{\mathcal{R}}\rangle +\displaystyle \frac{\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle }{2{\mathcal{R}}}[\displaystyle \frac{1-2{\mathcal{R}}}{1-{\mathcal{R}}}+\displaystyle \frac{2}{1+{({\mathcal{R}}/{{\mathcal{R}}}_{2})}^{n}}].\end{eqnarray} \tag{ 69 }$$

In the limit ${\mathcal{R}}\ll \{{{\mathcal{R}}}_{1},{{\mathcal{R}}}_{2}\}$, these expressions become

$$\begin{eqnarray}{D}_{{\mathcal{R}}} & \to & 6A({\mathcal{E}})[{(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{2}-{(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{1}})}^{2}],\\ {D}_{{\mathcal{R}}{\mathcal{R}}} & \to & 2A({\mathcal{E}}){\mathcal{R}}{(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{2}.\end{eqnarray} \tag{ 70 }$$

The ${\mathcal{R}}$-directed flux is

$$\begin{eqnarray}&&{\varphi }_{{\mathcal{R}}}=-{D}_{{\mathcal{R}}{\mathcal{R}}}\displaystyle \frac{\partial f}{\partial {\mathcal{R}}}-{D}_{{\mathcal{R}}}f\end{eqnarray} \tag{ 71 }$$

which, in the small-${\mathcal{R}}$ limit, becomes

$$\begin{eqnarray}{\varphi }_{{\mathcal{R}}} & \to & -2A({\mathcal{E}}){(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{2}(3\lambda f+{\mathcal{R}}\displaystyle \frac{\partial f}{\partial {\mathcal{R}}}),\\ \lambda & \equiv & 1-{(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{{\mathcal{R}}}_{1}})}^{2}.\end{eqnarray} \tag{ 72 }$$

Steady-state solutions can be characterized by either a constant, or a zero, flux. Setting ${\varphi }_{{\mathcal{R}}}=0$ yields

$$\begin{eqnarray}&&f({\mathcal{R}})=f({{\mathcal{R}}}_{2}){(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{-3\lambda },\ {\mathcal{R}}\ll \{{{\mathcal{R}}}_{1},{{\mathcal{R}}}_{2}\}.\end{eqnarray} \tag{ 73 }$$

When ${{\mathcal{R}}}_{1}\lt {{\mathcal{R}}}_{2}$, $\lambda \lt 0$ and the solution decreases toward ${\mathcal{R}}=0$; the reverse is true when ${{\mathcal{R}}}_{1}\gt {{\mathcal{R}}}_{2}$. Setting ${{\mathcal{R}}}_{1}={{\mathcal{R}}}_{2}$ yields $\lambda =0$ and $f({\mathcal{R}})=\mathrm{const}.$, the "zero-drift" solution.

A steady-state solution with constant but nonzero flux has the small-${\mathcal{R}}$ form

$$\begin{eqnarray}&&f({\mathcal{R}})={(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{-2}[\displaystyle \frac{1}{2-3\lambda }\displaystyle \frac{{\varphi }_{{\mathcal{R}}}}{2A}+{C}_{\lambda }{(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{2-3\lambda }]\end{eqnarray} \tag{ 74 }$$

for $\lambda \ne 2/3$, with ${C}_{\lambda }$ an integration constant. For $\lambda =2/3$,

$$\begin{eqnarray}&&f({\mathcal{R}})={(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{-2}[-\displaystyle \frac{{\varphi }_{{\mathcal{R}}}}{2A}\mathrm{log}(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})+{C}_{\frac{2}{3}}].\end{eqnarray} \tag{ 75 }$$

We compute ${C}_{\lambda }$ by requiring f to fall to zero at ${\mathcal{R}}\equiv {{\mathcal{R}}}_{0}$ ("empty loss cone"). The results are

$$\begin{eqnarray}&&f({\mathcal{R}})=f({{\mathcal{R}}}_{2}){(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{-2}\displaystyle \frac{{({\mathcal{R}}/{{\mathcal{R}}}_{0})}^{q}-1}{{({{\mathcal{R}}}_{2}/{{\mathcal{R}}}_{0})}^{q}-1},\ \ \ q\equiv 2-3\lambda ,\end{eqnarray} \tag{ 76a }$$

$$\begin{eqnarray}&&{\varphi }_{{\mathcal{R}}}({\mathcal{E}})=-2{qA}({\mathcal{E}})\displaystyle \frac{f({{\mathcal{R}}}_{2})}{1-{({{\mathcal{R}}}_{2}/{{\mathcal{R}}}_{0})}^{q}}\end{eqnarray} \tag{ 76b }$$

for $\lambda \ne 2/3$, and

$$\begin{eqnarray}&&f({\mathcal{R}})=f({{\mathcal{R}}}_{2}){(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{-2}\displaystyle \frac{\mathrm{log}({\mathcal{R}}/{{\mathcal{R}}}_{0})}{\mathrm{log}({{\mathcal{R}}}_{2}/{{\mathcal{R}}}_{0})},\end{eqnarray} \tag{ 77a }$$

$$\begin{eqnarray}&&{\varphi }_{{\mathcal{R}}}({\mathcal{E}})=-2A({\mathcal{E}})\displaystyle \frac{f({{\mathcal{R}}}_{2})}{\mathrm{log}({{\mathcal{R}}}_{2}/{{\mathcal{R}}}_{0})}\end{eqnarray} \tag{ 77b }$$

for $\lambda =2/3$.

At most energies, ${{\mathcal{R}}}_{0}\ll {{\mathcal{R}}}_{2}$. Assuming this inequality, f and ${\varphi }_{{\mathcal{R}}}$ have the following forms, depending on the value of λ:

• 1.  
  $\lambda \gt 2/3$, i.e., ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}\gt \sqrt{3}$. Defining $p\equiv 3\lambda -2\gt 0$,

  $$\begin{eqnarray}\displaystyle \frac{f({\mathcal{R}})}{f({{\mathcal{R}}}_{2})} & \approx & {(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{\mathcal{R}}})}^{2}[1-{(\displaystyle \frac{{{\mathcal{R}}}_{0}}{{\mathcal{R}}})}^{p}],\\ {\varphi }_{{\mathcal{R}}} & \approx & -2pA({\mathcal{E}})f({{\mathcal{R}}}_{2}).\end{eqnarray} \tag{ 78 }$$
• 2.  
  $\lambda \lt 2/3$, i.e., ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}\lt \sqrt{3}$. In terms of $q\equiv 2-3\lambda \gt 0$,

  $$\begin{eqnarray}\displaystyle \frac{f({\mathcal{R}})}{f({{\mathcal{R}}}_{2})} & \approx & {(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{\mathcal{R}}})}^{2}[{(\displaystyle \frac{{\mathcal{R}}}{{{\mathcal{R}}}_{2}})}^{q}-{(\displaystyle \frac{{{\mathcal{R}}}_{0}}{{{\mathcal{R}}}_{2}})}^{q}],\\ {\varphi }_{{\mathcal{R}}} & \approx & -2qA({\mathcal{E}}){(\displaystyle \frac{{{\mathcal{R}}}_{0}}{{{\mathcal{R}}}_{2}})}^{q}f({{\mathcal{R}}}_{2}).\end{eqnarray} \tag{ 79 }$$

  The "zero-drift" case has ${{\mathcal{R}}}_{1}={{\mathcal{R}}}_{2}$, $\lambda =0$, q = 2.
• 3.  
  $\lambda =2/3$, i.e., ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}=\sqrt{3}$ :

  $$\begin{eqnarray}\displaystyle \frac{f({\mathcal{R}})}{f({{\mathcal{R}}}_{2})} & \approx & {(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{\mathcal{R}}})}^{2}\displaystyle \frac{\mathrm{log}({\mathcal{R}}/{{\mathcal{R}}}_{0})}{\mathrm{log}({{\mathcal{R}}}_{2}/{{\mathcal{R}}}_{0})},\\ {\varphi }_{{\mathcal{R}}} & \approx & -\displaystyle \frac{2A({\mathcal{E}})f({{\mathcal{R}}}_{2})}{\mathrm{log}({{\mathcal{R}}}_{2}/{{\mathcal{R}}}_{0})}.\end{eqnarray} \tag{ 80 }$$

It is interesting to compare the expressions for the flux to those that would obtain in the absence of anomalous relaxation. Repeating the analysis with ${g}_{1}={g}_{2}=1$, we find:

$$\begin{eqnarray}{D}_{{\mathcal{R}}} & = & 0,\ \ {D}_{{\mathcal{R}}{\mathcal{R}}}=2A({\mathcal{E}}){\mathcal{R}}(1-{\mathcal{R}}),\\ {\varphi }_{{\mathcal{R}}} & = & -2A({\mathcal{E}}){\mathcal{R}}(1-{\mathcal{R}})\displaystyle \frac{\partial f}{\partial {\mathcal{R}}},\end{eqnarray} \tag{ 81 }$$

so that the steady state is characterized by

$$\begin{eqnarray}&&2A({\mathcal{E}}){\mathcal{R}}(1-{\mathcal{R}})\displaystyle \frac{\partial f}{\partial {\mathcal{R}}}+{\varphi }_{{\mathcal{R}}}=0.\end{eqnarray} \tag{ 82 }$$

The solution is

$$\begin{eqnarray}&&\displaystyle \frac{f({\mathcal{R}})}{f({{\mathcal{R}}}_{2})}=\displaystyle \frac{\mathrm{log}(\displaystyle \frac{{\mathcal{R}}}{1-{\mathcal{R}}})-\mathrm{log}(\displaystyle \frac{{{\mathcal{R}}}_{0}}{1-{{\mathcal{R}}}_{0}})}{\mathrm{log}(\displaystyle \frac{{{\mathcal{R}}}_{2}}{1-{\mathcal{R}}})-\mathrm{log}(\displaystyle \frac{{{\mathcal{R}}}_{0}}{1-{{\mathcal{R}}}_{0}})}\end{eqnarray} \tag{ 83 }$$

with flux

$$\begin{eqnarray}&&{\varphi }_{{\mathcal{R}}}=-2A({\mathcal{E}})\displaystyle \frac{f({{\mathcal{R}}}_{2})}{\mathrm{log}(\displaystyle \frac{{{\mathcal{R}}}_{2}}{1-{{\mathcal{R}}}_{2}})-\mathrm{log}(\displaystyle \frac{{{\mathcal{R}}}_{0}}{1-{{\mathcal{R}}}_{0}}).}\end{eqnarray} \tag{ 84 }$$

Again supposing that ${{\mathcal{R}}}_{0}\ll {{\mathcal{R}}}_{2}$, these expressions become:

$$\begin{eqnarray}\displaystyle \frac{f({\mathcal{R}})}{f({{\mathcal{R}}}_{2})} & \approx & \displaystyle \frac{\mathrm{log}{\mathcal{R}}-\mathrm{log}{{\mathcal{R}}}_{0}}{\mathrm{log}{{\mathcal{R}}}_{2}-\mathrm{log}{{\mathcal{R}}}_{0}},\\ {\varphi }_{{\mathcal{R}}} & \approx & -\displaystyle \frac{2A({\mathcal{E}})f({{\mathcal{R}}}_{2})}{\mathrm{log}({{\mathcal{R}}}_{2}/{{\mathcal{R}}}_{0})}.\end{eqnarray} \tag{ 85 }$$

Define a "reduction factor," η, as the ratio of this flux to the flux that would obtain in the presence of anomalous relaxation, assuming the same value for $f({{\mathcal{R}}}_{2})$. Then

$$\begin{eqnarray*}\eta \approx \{\begin{array}{lc}p\mathrm{log}({{\mathcal{R}}}_{2}/{{\mathcal{R}}}_{0}) & {\rm{if}}\;{\rm{}}{{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}\gt \sqrt{3}\\ q{({{\mathcal{R}}}_{0}/{{\mathcal{R}}}_{2})}^{q}\mathrm{log}({{\mathcal{R}}}_{2}/{{\mathcal{R}}}_{0}) & {\rm{if}}\;{\rm{}}{{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}\lt \sqrt{3}\\ 1 & {\rm{if}}\;{\rm{}}{{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}=\sqrt{3}.\end{array}\end{eqnarray*}$$

Figures 10 and 11 plot more accurate expressions for $f({\mathcal{R}})$ and η, computed using Equations (68a) and (71), without assuming the smallness of ${\mathcal{R}}$ or ${{\mathcal{R}}}_{0}$.

The "zero-drift" solution has ${{\mathcal{R}}}_{1}={{\mathcal{R}}}_{2}$, $\lambda =0$ and q = 2. For these values,

$$\begin{eqnarray}&&\eta =2{(\displaystyle \frac{{{\mathcal{R}}}_{0}}{{{\mathcal{R}}}_{2}})}^{2}\mathrm{log}(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{{\mathcal{R}}}_{0}}).\end{eqnarray} \tag{ 86 }$$

For ${{\mathcal{R}}}_{0}/{{\mathcal{R}}}_{2}=\{0.5,0.1,0.01,0.001\}$, $\eta \approx \{0.35,0.046,9.2\times {10}^{-4},1.4\times {10}^{-5}\}$.

The value of ${{\mathcal{R}}}_{1}/{{\mathcal{R}}}_{2}$ favored in the numerical experiments was ∼0.8, implying $\lambda \approx -0.56$ and $q\approx 3.7$. For these values,

$$\begin{eqnarray}&&\eta \approx 3.7{(\displaystyle \frac{{{\mathcal{R}}}_{0}}{{{\mathcal{R}}}_{2}})}^{3.7}\mathrm{log}(\displaystyle \frac{{{\mathcal{R}}}_{2}}{{{\mathcal{R}}}_{0}}).\end{eqnarray} \tag{ 87 }$$

For ${{\mathcal{R}}}_{0}/{{\mathcal{R}}}_{2}=\{0.5,0.1,0.01,0.001\}$, $\eta \approx \{0.20,1.7\times {10}^{-3},6.8\times {10}^{-7},2.0\times {10}^{-10}\}$.

Next we identify $\{{w}_{1},{w}_{2}\}$ with $\{{h}_{1},{h}_{2}\}$ given by Equations (42) and (41):

$$\begin{eqnarray}{h}_{1}({\mathcal{E}},{\mathcal{R}}) & = & [1+\displaystyle \frac{2(1-{\mathcal{R}})}{1-2{\mathcal{R}}}{(\displaystyle \frac{{{\mathcal{R}}}_{3}}{{\mathcal{R}}})}^{2}]\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{3}^{2}}{{{\mathcal{R}}}^{2}}),\\ {h}_{2}({\mathcal{E}},{\mathcal{R}}) & = & \mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}}).\end{eqnarray} \tag{ 88 }$$

The flux coefficients are

$$\begin{eqnarray}{D}_{{\mathcal{R}}} & = & 2A({\mathcal{E}})[1-2{\mathcal{R}}+2(1-{\mathcal{R}}){(\displaystyle \frac{{{\mathcal{R}}}_{4}}{{\mathcal{R}}})}^{2}]\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}})\\ & & -2A({\mathcal{E}})[1-2{\mathcal{R}}+2(1-{\mathcal{R}}){(\displaystyle \frac{{{\mathcal{R}}}_{3}}{{\mathcal{R}}})}^{2}]\\ & & \times \mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{3}^{2}}{{{\mathcal{R}}}^{2}}),\end{eqnarray} \tag{ 89a }$$

$$\begin{eqnarray}&&{D}_{{\mathcal{R}}{\mathcal{R}}}=2A({\mathcal{E}}){\mathcal{R}}(1-{\mathcal{R}})\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}})\end{eqnarray} \tag{ 89b }$$

and

$$\begin{eqnarray}&&{D}_{{\mathcal{R}}}=-\langle {\rm{\Delta }}{\mathcal{R}}\rangle +\displaystyle \frac{\langle {({\rm{\Delta }}{\mathcal{R}})}^{2}\rangle }{2{\mathcal{R}}}(\displaystyle \frac{1-2{\mathcal{R}}}{1-{\mathcal{R}}}+2\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}}).\end{eqnarray} \tag{ 90 }$$

In the limit ${\mathcal{R}}\ll \{{{\mathcal{R}}}_{3},{{\mathcal{R}}}_{4}\}$, these expressions become

$$\begin{eqnarray}{D}_{{\mathcal{R}}} & \to & 4A({\mathcal{E}}){(\displaystyle \frac{{{\mathcal{R}}}_{4}}{{\mathcal{R}}})}^{2}\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}})\\ & & \times [1-\displaystyle \frac{{{\mathcal{R}}}_{3}^{2}}{{{\mathcal{R}}}_{4}^{2}}\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{3}^{2}-{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}})],\\ {D}_{{\mathcal{R}}{\mathcal{R}}} & \to & 2A({\mathcal{E}}){\mathcal{R}}\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}}).\end{eqnarray} \tag{ 91 }$$

In the same limit, the ${\mathcal{R}}$-directed flux is

$$\begin{eqnarray}{\varphi }_{{\mathcal{R}}} & \to & -2A({\mathcal{E}})\mathrm{exp}(-\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}})\\ & & \times \{2{(\displaystyle \frac{{{\mathcal{R}}}_{4}}{{\mathcal{R}}})}^{2}[1-\displaystyle \frac{{e}^{-\delta {{\mathcal{R}}}_{3}^{2}/{{\mathcal{R}}}^{2}}}{1-\delta }]f+{\mathcal{R}}\displaystyle \frac{\partial f}{\partial {\mathcal{R}}}\},\\ \delta & \equiv & 1-{(\displaystyle \frac{{{\mathcal{R}}}_{4}}{{{\mathcal{R}}}_{3}})}^{2}.\end{eqnarray} \tag{ 92 }$$

Setting ${\varphi }_{{\mathcal{R}}}=0$ yields

$$\begin{eqnarray}\displaystyle \frac{\mathrm{log}f({\mathcal{R}})}{\mathrm{log}f({{\mathcal{R}}}_{4})} & = & \displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}}-1+\displaystyle \frac{{{\mathcal{R}}}_{3}^{2}}{{{\mathcal{R}}}_{3}^{2}-{{\mathcal{R}}}_{4}^{2}}[\mathrm{exp}(\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}-{{\mathcal{R}}}_{3}^{2}}{{{\mathcal{R}}}^{2}})\\ & & -\mathrm{exp}(\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}-{{\mathcal{R}}}_{3}^{2}}{{{\mathcal{R}}}_{4}^{2}})],\ {\mathcal{R}}\ll \{{{\mathcal{R}}}_{3},{{\mathcal{R}}}_{4}\}.\end{eqnarray} \tag{ 93 }$$

For ${{\mathcal{R}}}_{3}\lt {{\mathcal{R}}}_{4}$, the dominant terms imply

$$\begin{eqnarray}&&f({\mathcal{R}})\to f({{\mathcal{R}}}_{4})\mathrm{exp}[{\delta }^{-1}\mathrm{exp}(\alpha \displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}})]\end{eqnarray} \tag{ 94 }$$

where $\delta =1-{{\mathcal{R}}}_{4}^{2}/{{\mathcal{R}}}_{3}^{3}\lt 0$, $\alpha =1-{{\mathcal{R}}}_{3}^{2}/{{\mathcal{R}}}_{4}^{2}\gt 0$, hence f drops very rapidly to zero below ${\mathcal{R}}={{\mathcal{R}}}_{4}$. Whereas for ${{\mathcal{R}}}_{3}\gt {{\mathcal{R}}}_{4}$,

$$\begin{eqnarray}&&f({\mathcal{R}})\to f({{\mathcal{R}}}_{4}){e}^{{{\mathcal{R}}}_{4}^{2}/{{\mathcal{R}}}^{2}},\end{eqnarray} \tag{ 95 }$$

a rapid rise toward ${\mathcal{R}}=0$. This is qualitatively the same behavior as in the power-law case when ${\varphi }_{{\mathcal{R}}}=0$.

In the case of a constant but nonzero flux, only the "zero-drift" solution (with ${{\mathcal{R}}}_{3}={{\mathcal{R}}}_{4}$) can be expressed in terms of simple functions:

$$\begin{eqnarray}f({\mathcal{R}}) & \approx & -\displaystyle \frac{{\varphi }_{{\mathcal{R}}}}{4A}[{E}_{i}(\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}_{0}^{2}})-{E}_{i}(\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}^{2}})],\\ & & {\mathcal{R}}\ll {{\mathcal{R}}}_{3},{{\mathcal{R}}}_{4}\end{eqnarray} \tag{ 96a }$$

$$\begin{eqnarray}&&\displaystyle \frac{f({\mathcal{R}})}{f({{\mathcal{R}}}_{4})}\approx \displaystyle \frac{{E}_{i}({{\mathcal{R}}}_{4}^{2}/{{\mathcal{R}}}_{0}^{2})-{E}_{i}({{\mathcal{R}}}_{4}^{2}/{{\mathcal{R}}}^{2})}{{E}_{i}({{\mathcal{R}}}_{4}^{2}/{{\mathcal{R}}}_{0}^{2})-{E}_{i}(1)},\end{eqnarray} \tag{ 96b }$$

$$\begin{eqnarray}&&-\displaystyle \frac{{\varphi }_{{\mathcal{R}}}}{4A}\approx f({{\mathcal{R}}}_{4}){[{E}_{i}(\displaystyle \frac{{{\mathcal{R}}}_{4}^{2}}{{{\mathcal{R}}}_{0}^{2}})-{E}_{i}(1)]}^{-1}.\end{eqnarray} \tag{ 96c }$$

In these expressions, E_i is the exponential function. The reduction factor defined above becomes, in this case,

$$\begin{eqnarray}&&\eta =2\displaystyle \frac{{{\mathcal{R}}}_{4}}{{{\mathcal{R}}}_{0}}\mathrm{log}(\displaystyle \frac{{{\mathcal{R}}}_{4}}{{{\mathcal{R}}}_{0}}){e}^{-{{\mathcal{R}}}_{4}/{{\mathcal{R}}}_{0}}\end{eqnarray} \tag{ 97 }$$

For ${{\mathcal{R}}}_{0}/{{\mathcal{R}}}_{4}=\{0.5,0.1,0.01\}$, $\eta \approx \{0.38,2.1\times {10}^{-3},3.4\times {10}^{-41}\}$.