The statistical mechanics of relativistic orbits around a massive black hole

TL;DR

This paper develops a statistical mechanics framework to analyze how correlated Gaussian noise influences relativistic stellar orbits around massive black holes, revealing a barrier effect in angular momentum evolution due to Schwarzschild precession.

Contribution

It introduces a phase-averaged stochastic Hamiltonian approach and derives the Fokker-Planck equation for correlated Gaussian noise, highlighting the impact of noise smoothness on orbital dynamics.

Findings

Angular momentum evolution is suppressed below a critical value j_b.

A Schwarzschild precession-induced barrier in angular momentum is identified.

The barrier's significance depends on the noise's temporal smoothness and frequency.

Abstract

Stars around a massive black hole (MBH) move on nearly fixed Keplerian orbits, in a centrally-dominated potential. The random fluctuations of the discrete stellar background cause small potential perturbations, which accelerate the evolution of orbital angular momentum by resonant relaxation. This drives many phenomena near MBHs, such as extreme mass-ratio gravitational wave inspirals, the warping of accretion disks, and the formation of exotic stellar populations. We present here a formal statistical mechanics framework to analyze such systems, where the background potential is described as a correlated Gaussian noise. We derive the leading order, phase-averaged 3D stochastic Hamiltonian equations of motion, for evolving the orbital elements of a test star, and obtain the effective Fokker-Planck equation for a general correlated Gaussian noise, for evolving the stellar distribution function. We show that the evolution of angular momentum depends critically on the temporal smoothness of the background potential fluctuations. Smooth noise has a maximal variability frequency $\nu_{\max}$. We show that in the presence of such noise, the evolution of the normalized angular momentum $j=\sqrt{1-e^{2}}$ of a relativistic test star, undergoing Schwarzschild (in-plane) General Relativistic precession with frequency $\nu_{GR}/j^{2}$, is exponentially suppressed for $j<j_{b}$, where $\nu_{GR}/j_{b}^{2}\sim\nu_{\max}$, due to the adiabatic invariance of the precession against the slowly varying random background torques. This results in an effective Schwarzschild precession-induced barrier in angular momentum. When $j_{b}$ is large enough, this barrier can have significant dynamical implications for processes near the MBH.

Equations

$$H\left(a,J,J_z,\varphi ,\psi ,t\right)=H_K\left(a\right)+H_{GR}\left(a,J\right)+U\left(a,J,J_z,\varphi ,\psi ,t\right),$$

$$U\left(a,J,J_z,\varphi ,\psi ,t\right)=\sum _{l=0}^{\infty }{H}_l\left(a,J,J_z,\varphi ,\psi ,t\right),$$

$$H_l\left(a,J,J_z,\varphi ,\psi ,t\right)$$

$$h_{nm}^l\left(a,J,t\right)$$

$$⟨h_{00}^0\left(a,J,t\right)⟩=4\pi N{\tau }_0\int n\left(a_k,e_k\right)d{a}_{k}d{e}_k{⟨\frac{a}{r}{K}_0\left(r_k/r\right)⟩}_{\circlearrowright },$$

$$⟨h_{nm}^l\left(a,J,t\right)⟩=0,l>0,$$

$$⟨{\eta }_{nm}^l\left(a,J,t\right){\eta }_{n^{\prime }{m}^{\prime }}^{l^{\prime }\ast }\left(a^{\prime },J^{\prime },t^{\prime }\right)⟩$$

$$H\left(a,J,J_z,\varphi ,t\right)$$

$${\eta }_{\pm 11}^1\left(a,J,t\right)$$

$${\eta }_{\pm 1-1}^1\left(a,J,t\right)$$

$${\eta }_{\pm 10}^1\left(a,J,t\right)$$

$$⟨{\eta }_i\left(t\right){\eta }_j\left(t^{\prime }\right)⟩={\delta }_{ij}C\left(t-t^{\prime }\right),$$

$$T_{\mathrm{coh}}={\int }_0^{\infty }C\left(t\right)dt.$$

$$H_{\eta }=-{\nu }_{GR}{J}_c^2/J+⟨h_{00}^0\left(J,a,t\right)⟩+{\eta }_{00}^0\left(J,a,t\right)+\sqrt{2C^1\left(a,J\right)}{\hat{e}}_{\psi }(\mathbf{\Pi })\cdot \mathbf{\eta }\left(t\right),$$

$${\hat{e}}_{\psi }(\mathbf{\Pi })\equiv \sin \psi {\hat{e}}_{\varphi }(\varphi )-\cos \psi {\hat{e}}_u(\varphi ,u).$$

$$\hat{e}_{J}=\left(\begin{array}[]{c}\sqrt{1-u^{2}}\cos\phi\\ \sqrt{1-u^{2}}\sin\phi\\ u\end{array}\right);\;\hat{e}_{\phi}=\left(\begin{array}[]{c}-\sin\phi\\ \cos\phi\\ 0\end{array}\right);\;\hat{e}_{u}=\left(\begin{array}[]{c}-u\cos\phi\\ -u\sin\phi\\ \sqrt{1-u^{2}}\end{array}\right),$$

$$H_{\eta }^{GR}=-{\nu }_{GR}{J}_c/j+{\tau }_N^1(j){\hat{e}}_{\psi }(\mathbf{\Pi })\cdot \mathbf{\eta }\left(t\right),$$

$$\dot{\mathbf{J}}$$

$$\dot{\psi }$$

$$\dot{x}$$

$$\dot{\psi }$$

$${\mathbf{\nu }}_j\left(j,\mathbf{\Pi }\right)$$

$${\mathbf{\nu }}_{\varphi }\left(j,\mathbf{\Pi }\right)$$

$${\mathbf{\nu }}_u\left(j,\mathbf{\Pi }\right)$$

$${\mathbf{\nu }}_{\psi }\left(j,\mathbf{\Pi }\right)$$

$$\frac{\partial }{\partial t}\phi \left(j,t\right)$$

$$\frac{\partial }{\partial t}P\left(j,t\right)=-\frac{\partial }{\partial j}⟨{\mathbf{\nu }}_j\left(j,\mathbf{\Pi }\left(t\right)\right)\cdot \mathbf{\eta }\left(t\right)\phi \left(j,t\right)⟩.$$

$$⟨{\eta }_i\left(t\right)R\left[\mathbf{\eta }(t^{\prime })\right]⟩=\int ds⟨{\eta }_i\left(t\right){\eta }_j\left(s\right)⟩⟨\frac{\delta R\left[\mathbf{\eta }(t^{\prime })\right]}{\delta {\eta }_j\left(s\right)}⟩.$$

$$\frac{\partial }{\partial t}P\left(j,t\right)$$

$${\nabla }_{\eta }\left(t\right)\equiv \left(\frac{\delta }{\delta {\eta }_1\left(t\right)},\frac{\delta }{\delta {\eta }_2\left(t\right)},\frac{\delta }{\delta {\eta }_3\left(t\right)}\right).$$

$$\frac{\partial }{\partial t}P\left(j,t\right)$$

$$x\left(t\right)$$

$$\psi \left(t\right)$$

$${\nabla }_{\eta }\left(t^{\prime }\right)x\left(t\right)$$

Full text

††: Class. Quantum Grav.

Keywords: Black holes, Galactic nuclei, Stellar dynamics and kinematics, perturbation theory

1 Introduction

Stars around a massive black hole (MBH) move on nearly fixed Keplerian orbits, in a centrally-dominated potential. The discrete stellar background and its random fluctuations in time, add only small perturbations to the MBH potential, which persist on a coherence timescale $T_{\mathrm{coh}}$. However, these perturbations induce a residual random torque of magnitude $\tau_{N}\propto\sqrt{N}$, where $N$ is the number of stars on the relevant scale, which leads to rapid evolution of the orbital angular momentum by the process of resonant relaxation [1]. The torques persist coherently until the background is randomized by the slow drift away from the fixed Keplerian orbits, due for example to the mean potential of the distributed mass (mass precession), General Relativistic (GR) precession, or ultimately by their response to the resonant torques themselves. Resonant relaxation is the driving force behind many interesting physical phenomena, which include the emission of gravitational waves from compact objects that spiral into the MBH [2], the warping of circum-nuclear gaseous or stellar disks [3, 4], and the orbital evolution of tidally captured stars [5, 6, 7]. Resonant relaxation can also be a major obstacle for attempts to test GR by observation of stars very near MBHs [8, 9].

The original treatment of resonant relaxation [1] assumed that a typical test orbit is similar to a typical background orbit, and therefore there is only one relevant time scale in the system—the coherence time $T_{\mathrm{coh}}$ of the background. However, eccentric ($e\to 1$) relativistic orbits undergo rapid Schwarzschild (in-plane) precession on a short timescale $(T_{GR}\propto j^{2}a^{5/2})<T_{\mathrm{coh}}(a)$, where $a$ is the orbital semi-major axis (sma) and $j=\sqrt{1-e^{2}}$ is the normalized orbital angular momentum. This deterministic orbital evolution changes the torques acting on the test star over one precession period. This motivated previous studies [1, 10] to adopt $T_{GR}$ as the effective stochastic coherence timescale of the random torquing process. These early studies noted that fast torquing by resonant relaxation is expected to become inefficient for eccentric orbits where $T_{GR}$ is much shorter than the coherence time of the background, and where gravitational wave emission is already dynamically significant. This would then allow inspiral by gravitational wave emission to proceed unimpeded. However, this analysis failed to take into account the periodic (i.e. regular rather than stochastic) nature of the precession. In fact, depending on the statistical properties of the background fluctuations, the effect of precession can go well beyond merely shortening the coherence time and suppressing resonant relaxation—it can actually lead to adiabatic invariance of the angular momentum. Indeed, direct $N$-body simulations [11, 12] revealed that the stochastic evolution of the angular momentum is restricted by relativistic precession at the so-called “Schwarzschild barrier“ [11]. However, a rigorous theoretical framework for describing these relativistic kinematics, in terms of kinetic and stochastic theory, is still lacking [13].

Our goal is to study relativistic dynamics near a MBH. We focus on the interplay between the deterministic relativistic precession and the stochastic fluctuations of the background potential. We derive here stochastic equations of motion (EOM), for evolving the orbital elements of a test star, and obtain the corresponding Fokker-Planck (FP) equation, for evolving the stellar distribution function.

In Section 2 we expand the orbit-averaged Hamiltonian of a test star, where the effect of the background potential is represented as a random correlated Gaussian noise. We then derive the stochastic EOM from the leading-order stochastic Hamiltonian. The effective FP equation for a general correlated Gaussian noise is then derived from the EOM in Section 3. In Section 4 we validate the derived FP equation by comparing its results to direct integration of the stochastic EOM. We then show that the evolution of the angular momentum depends not only on the values of the coherence time $T_{\mathrm{coh}}$ and the magnitude of the random torque $\tau_{N}$, as expected, but also critically on the smoothness (differentiability) of the background noise. Specifically, since a smooth background noise has a finite maximal variability frequency $\nu_{\max}$, the evolution of the test star’s angular momentum is exponentially suppressed for $j<j_{b}(a)$, where $T_{GR}(a,j_{b})=2\pi/\nu_{\max}(a)$. This results in an effective barrier in angular momentum due to adiabatic invariance induced by the relativistic precession. We argue that when $j_{b}$ is large enough to be relevant for stable relativistic orbits, this constitutes a Schwarzschild precession-induced barrier. We discuss the implications and summarize in section 5.

2 The stochastic Hamiltonian

Let $M_{\bullet}$ be the mass of the massive black hole (MBH), $M$ the mass of the test star, and $\nu_{r}=\sqrt{GM_{\bullet}/a^{3}}$ its Keplerian orbital frequency around the MBH, neglecting the mass of the background stars. The star’s orbital angular momentum is $J=J_{c}\sqrt{1-e^{2}}$, where $e$ is the eccentricity and $J_{c}=\sqrt{GM_{\bullet}a}$ is the maximal (circular) specific angular momentum for a given sma, $a$, i.e. the star’s specific binding energy in this limit is $E=GM_{\bullet}/2a$. The Euler angles $\boldsymbol{\Pi}=\left(\psi,\theta,\phi\right)$ describing the orbit’s orientation in a fixed inertial reference frame are linearly related to the standard orbital elements: the longitude of the ascending node is $\Omega=\phi+\pi/2$, the inclination angle is $i=\theta$ and the argument of pericenter is $\omega=\psi-\pi/2$. Note that neither the Euler angles nor the orbital elements are the canonical coordinates for the Hamiltonian.

The secular (orbit-averaged) Hamiltonian for a single star, in terms of canonical action-angle coordinates [14, 15, 16] $[(J,\psi),\,(J_{z},\phi),\,(I,w)]$, where, $J_{z}=J\cos\theta$ is the angular momentum in the $z$ direction and $I=\sqrt{GM_{\bullet}a}$, is

[TABLE]

where due to the orbit-averaging over the mean anomaly $w$, the action $I$, and therefore $a$, are constant. $H_{K}=-J_{c}\nu_{r}\left(a\right)/2$ is the Keplerian term, $H_{GR}=-\nu_{GR}J_{c}^{2}/J$ with $\nu_{GR}=3\nu_{r}\left(a\right)GM_{\bullet}/\left(c^{2}a\right)$ is the relativistic 1st post-Newtonian (PN) term, and $U\left(a,J,J_{z},\phi,\psi,t\right)$ is the orbit-averaged potential due to the background stars.

2.1 The potential of background stars

The time-dependent, randomly fluctuating stellar background potential can be expanded in Wigner rotation matrices $D_{mn}^{l}\left(\boldsymbol{\Pi}\right)$ (Bar-Or & Alexander 2014, in prep.). These are used here to describe the 3D position of an orbit, as given by its orbital elements, thereby allowing to generalize the standard Legendre expansion of the Newtonian potential of a collection of point masses, to that of objects that have an additional degree of freedom. This allows the description of both the orientation of the orbital planes, and the orientation of the extended orbit-averaged Keplerian ellipses in them111A common use of the Wigner matrices is to describe the rotation of point particles, where the additional degree of freedom is their intrinsic spin.,

[TABLE]

where the $l$’th degree spherical harmonic is further decomposed to $(2l+1)^{2}$ multipoles

[TABLE]

and where the terms

[TABLE]

depend on the distribution of the stellar background, which is a function of time through the evolution of the orbital elements of the $N$ stars, and is represented below by a noise model. Note that for each $(l,n,m)$ term, the time-dependence of the stellar background enters only through $h_{nm}^{l}$. The scale of the torque $\tau_{0}\left(a\right)=J_{c}\nu_{r}Q^{-1}$ due to a single background star is inversely proportional to the MBH-star mass ratio $Q=M_{\bullet}/M$ (a single-mass population assumed). The usual $\left(r_{<}/r_{>}\right)^{l}/r_{>}$ min-max terms that appear in the expansion of the potential to Legendre polynomials are rewritten in terms of the dimensionless functions $K_{l}\left(x\right)\equiv\min\left(x,1\right)^{2l+1}/x^{l+1}.$ The average $\left\langle\dots\right\rangle_{\circlearrowright}$ denotes orbit-averaging over the mean anomalies of both the test star and the $k$’th background star.

2.2 The background potential in the stochastic limit

A key assumption is that correlations between any two stars in the system are short-lived, and have a negligible affect on the long-term evolution of the system. Since each $h_{nm}^{l}\left(a,J,t\right)$ term results from the superposed gravitational forces by $N\gg 1$ background stars, we invoke the central limit theorem and assume that it can be described by time-dependent Gaussian random variables $\eta_{nm}^{l}\equiv h_{nm}^{l}-\left\langle h_{nm}^{l}\right\rangle$, which have zero mean and are therefore completely described by their 2-point correlation functions (Eq. 7). Therefore, although the system is completely deterministic, the time-dependent terms in the Hamiltonian can be considered as a time-correlated background “noise”. The derived EOM can then be interpreted as a set of non-linear Langevin equations (e.g. [17, Ch. 3]).

We assume that the stellar background is on average isotropic and stationary, so that only the spherical monopole component ($l=0$) has a non-zero mean,

[TABLE]

[TABLE]

where $\left\langle\ldots\right\rangle$ denotes the mean over all initial conditions and realizations of the background, and $n(a,e)$ is the number density of stars in $(a,e)$ phase space. The 2-point correlation functions $C_{n,n^{\prime}}^{l}$ are then simplified by the orthonormality of the Wigner $D$-matrices (Eq. 4)

[TABLE]

This expansion of the Hamiltonian into $(l,n,m)$ terms has the useful property of decoupling the angular dependence of the test orbit from that of the background for each term separately (Eq. 3). Note that since the auto-correlations $\left\langle\eta_{nm}^{l}\left(a,J,t\right)\eta_{nm}^{l^{\prime}*}\left(a,J,t\right)\right\rangle$ are independent of $m$, all orders of $m$ for a given $l$ and $n$ contribute equally to the Hamiltonian and should be taken into account.

2.3 First order ($l=1$) Hamiltonian

The $l\leq 1$ terms in the Hamiltonian are

[TABLE]

In general the $l=1$ term of the Hamiltonian has $\left(2l+1\right)^{2}=9$ real-valued $h_{nm}^{1}$ components, which correspond to 9 real-valued noise components222The $h_{00}^{0}$ term, which contributes only to the evolution of $\psi$ (precession due to the distributed stellar mass), is the only one with a non-zero mean. The mean $\left\langle h_{00}^{0}\right\rangle$ leads to a constant drift in $\psi$, while the noise associated with this term, $\eta_{00}^{0}=h_{00}^{0}-\left\langle h_{00}^{0}\right\rangle$ (i.e. fluctuations in the total stellar mass within a fixed radius), is a higher-order stochastic perturbation.. However, the symmetries of $h_{nm}^{1}$ imply that $\eta_{0m}^{1}=0$ and $\eta_{-1m}^{1}=\eta_{1m}^{1}$, leaving only three non-zero noise terms. We can therefore decompose the noise into three independent real Gaussian components $\boldsymbol{\tilde{\eta}}\left(a,J,t\right)=\left(\eta_{1}\left(a,J,t\right),\eta_{2}\left(a,J,t\right),\eta_{3}\left(a,J,t\right)\right)^{T}$

[TABLE]

Note that the noise components $\tilde{\eta}_{i}$ depend on the position of the test orbit in $(a,J)$ phase-space, and likewise, so do the $2$-point correlation functions. Here we simplify this general description by assuming that on the relevant timescales, which are shorter than the 2-body energy relaxation timescale, $a$ remains nearly constant, so the coherence scales can be evaluated at the fixed, initial $a$ value.

We further assume that at any given short time interval $(t,t+\mathrm{d}t)$ over which the test orbit’s angular momentum evolves over the interval $(J,J+\mathrm{d}J)$, the temporal correlations decay much faster than the angular momentum correlations, $\left\langle\eta\left(J,t\right)\eta\left(J,t+dt\right)\right\rangle\ll\left\langle\eta\left(J,t\right)\eta\left(J+dJ,t\right)\right\rangle$, where $dJ\sim\tau_{N}dt$. This assumption is motivated by the fact that $T_{\mathrm{coh}}<J_{c}/\tau_{N}$, whether the coherence time is determined by mass precession (e.g. [18]), or whether it is determined by the resonant torques themselves (as indicated by fully self-consistent numerical evolution simulations of circular orbits, Bar-Or & Alexander 2014, in prep.). This allows the separation $\boldsymbol{\tilde{\eta}}\left(a,J,t\right)\approx\sqrt{C^{1}\left(a,J\right)}\boldsymbol{\eta}\left(t\right)$, where the amplitude $C^{1}\left(a,J\right)=C_{11}^{1}\left(a,a,J,J,0\right)$ is normalized to reproduce the covariance matrix, and where $\boldsymbol{\eta}$ is a vector with three independent Gaussian components, with zero mean and an auto-correlation function (ACF)

[TABLE]

where $C(0)=1$. The ACF is characterized by its shape and magnitude. We define here the coherence time as

[TABLE]

Note that as defined, $T_{\mathrm{coh}}$ is half the total power of the noise.

The assumption of separability leads to a great simplification of the stochastic EOM, since the noise term $\boldsymbol{\eta}$ then enters as function of time only, without derivatives of the stochastic noise with respect to $J$. The $l\leq 1$ stochastic Hamiltonian then reads

[TABLE]

where

[TABLE]

and we introduce an orthonormal spherical coordinate system $\left(J,\phi,u\right)$, where $u\equiv\cos\theta$, with the associated unit vectors $\hat{e}_{i}=\left(\partial\mathbf{J}/\partial i\right)/\left|\partial\mathbf{J}/\partial i\right|$ for $i\in\left\{J,\phi,u\right\}$

[TABLE]

The noise term $\boldsymbol{\eta}$ can be viewed as a 3D vector in angular momentum space. In this sense, $\boldsymbol{\eta}$ plays a role similar to the “dipole” term in the Hamiltonian model of Merritt et al. [11] (see Eq. A6 there), which used the ansatz of an effective potential to represent the residual torque by the background stars. Here, in contrast, we derive the Hamiltonian from first principles without invoking an effective potential, and include the stochastic effect of the background stars as a formally defined noise term.

2.4 Integrating the stochastic equations of motion

We are interested here only in the coupling between the GR precession and the stochastic Newtonian perturbations of the background stars on the test orbit, we therefore omit the Newtonian spherical terms (mass precession terms), and $H_{\eta}$ reduces to

[TABLE]

where $\tau_{N}^{1}=\sqrt{2C^{1}\left(a,j\right)}\sim\sqrt{N(a)}\tau_{0}$ is the scale of the residual torque, which reflects the random contribution of the $N\left(a\right)$ stars on the relevant scale. Using static wires simulation similar to [19] but with a faster method [20], we calculated the total residual torque $\tau_{N}$ on a test wire in the $\hat{e}_{J}$ direction due to a cluster of $N$-wires. The total torque can be approximated by $\tau_{N}\left(a,j\right)=\sqrt{\left\langle\tau_{J}^{2}\right\rangle}\approx 0.28\sqrt{1-j}\sqrt{N\left(2a\right)}GM/a$. This is larger than the partial torque $\tau_{N}^{1}$ corresponding to the $l=1$ term, which rises asymptotically to $\approx 0.6\tau_{N}\left(a,j=0\right)$ as $j\to 0$ (by numeric integration of $C_{1,1}^{1}$, using Eq. 4). Nevertheless, we replace $\tau_{N}^{1}$ with $\tau_{N}$. This has the advantage of simplicity, and of reproducing the total power of the residual torques with the $l=1$ approximation, which may be a better model of the real physical system.

The EOM derived from $H_{\eta}^{GR}$ can be reduced to the compact form

[TABLE]

In this form the EOM are non-linear Langevin equations, which require a specified noise model $\boldsymbol{\eta}(t)$. The system can be integrated in time from specified initial conditions, for any given realization of the noise. The statistical properties of the system can then be obtained by repeating this for many noise realizations, in a Monte Carlo fashion. Alternatively, the evolution of the probability density function (PDF) of the test star in phase-space can be derived from the corresponding FP equation. This is a well-defined procedure in the Markovian limit, where $\boldsymbol{\eta}$ can be treated as an uncorrelated white noise on all physically relevant timescales. In the general non-Markovian case, it is not always possible to describe the dynamics by an effective FP equation. However, the system of interest here can be approximately described in this form (Section 3), as is validated by comparison to direct integration of the stochastic EOM (Section 4).

3 Effective diffusion with correlated noise

In the $N\gg 1$ limit, the fluctuating potential of the background stars can be treated as a stochastic noise. The EOM for $\phi$, $u=\cos\theta$, $\psi$ and total angular momentum $j=J/J_{c}$ are then derived from the stochastic Hamiltonian $H_{\eta}$ (Eq. 14),

[TABLE]

where the vectors $\boldsymbol{\nu}_{x}\left(j,\boldsymbol{\Pi}\right)$ express the torques in the $(j,\phi,u,\omega)$ coordinates, and transform the noise $\boldsymbol{\eta}$ to them from the $\mathbf{J}$ phase-space in which it is defined.

[TABLE]

with $\nu_{j}\left(j\right)=\left|\boldsymbol{\mathbf{\nu}}_{j}\left(j\right)\right|=\tau_{N}\left(J\right)/J_{c}$ and $\hat{e}_{\psi\psi}=\partial\hat{e}_{\psi}/\partial\psi$.

These EOM are a set of non-linear Langevin type equations. The corresponding FP equations, which describe the evolution of the PDF of the state variables, $P\left(j,\phi,u,\psi,t\right)$, can be derived in the Markovian limit (e.g. [17]). For simplicity, we focus here only on the PDF of the normalized angular momentum, $j$, marginalizing over all the other state variables.

For a given realization of the noise, the trajectory of a test star in phase-space, $(j(t),\boldsymbol{\Pi}(t)),$ is fully specified by Eqs. (20, 21) and the initial conditions. Therefore, its $j$-trajectory can be formally described by $\varphi\left(j,t\right)=\delta\left(j-j\left(t\right)\right)$, which obeys the continuity equation

[TABLE]

The PDF is then given by $P\left(j,t\right)=\left\langle\varphi\left(j,t\right)\right\rangle$, where $\left\langle\dots\right\rangle$ denotes the average over all realizations and all initial conditions of the noise. The evolution of $P\left(j,t\right)$ is determined by

[TABLE]

Since the noise term $\boldsymbol{\eta}\left(t\right)$ is a Gaussian with zero mean, Novikov’s theorem [21], which expresses the correlation of the noise with any functional $R\left[\boldsymbol{\eta}\right]$ of the noise, can be applied,

[TABLE]

Therefore, Eqs. (12) and (27) imply

[TABLE]

where we define for brevity the functional gradient operator

[TABLE]

By multiple application of the chain rule, Eq. (29) reads

[TABLE]

We integrate the EOM (Eqs. 20, 21) to obtain

[TABLE]

Thus, for $t^{\prime}\leq t$, the response functions for $x\in\left\{j,\phi,u,\psi\right\}$ are

[TABLE]

These functions describe the manner by which a small perturbation of the noise at an initial time $t^{\prime}$ is magnified at a later time $t$ by the accumulated differences in torques experienced along the perturbed and unperturbed trajectories, due to the fact that the noise amplitude ($\boldsymbol{\nu}_{x}$) depends on the state variables along the trajectory.

In general, Eq. (LABEL:eq:P_dot_Nov2) can be written in closed form (i.e. solely in terms of $P$) only when the response of the state variables to the noise (given by the response functions $\nabla_{\eta}\left(t^{\prime}\right)x\left(t\right)$) can be reduced to a form that involves only deterministic functions (i.e. a known function of $j(t),$ $t$ and $t^{\prime}$). Otherwise, the Novikov theorem has to be applied recursively, but closure is not guaranteed [22, 23].

To address this difficulty, we first consider the Markovian limit, where this problem is avoided. In this limit, the process depends only on the total power of the noise, and therefore can be presented by an appropriately normalized uncorrelated noise. Eq. (LABEL:eq:P_dot_Nov2) then reduces to the FP equation (section 3.1). We further show below, that in the non-Markovian case, which is relevant in the relativistic limit where the precession time is short, the response functions can be approximated by deterministic functions, which we then employ to convert Eq. (LABEL:eq:P_dot_Nov2) to an effective FP equation (section 3.2).

3.1 The Markovian Limit

In the limit where the correlation time is shorter than any other relevant timescale in the system, the noise can be regarded as an uncorrelated (white) noise with a correlation function $C\left(|t-t^{\prime}|\right)=2T_{\mathrm{coh}}\delta\left(t-t^{\prime}\right)$, and the process is Markovian333Note that in this singular case, $C(0)\to\infty$ is no longer normalized to $C(0)=1$ as previously assumed. This inconsistency is not a problem, because in this limit the dynamics depend only on the total power of the noise, $2T_{\mathrm{coh}}$ (Eq. 13), through the diffusion coefficient $D_{2}=2\tau_{N}^{2}T_{\mathrm{coh}}$.. In this limit, the response functions contribute to the integrals in Eq. (LABEL:eq:P_dot_Nov2) only at the limit $t^{\prime}=t$ and therefore

[TABLE]

where we used

[TABLE]

It then follows from

[TABLE]

that the FP equation is

[TABLE]

where $D_{2}(j)=D(\Delta j^{2})=2T_{\mathrm{coh}}\nu_{j}^{2}\left(j\right)$ is the usual 2nd order diffusion coefficient (DC). It is easy to verify that the zero-flux steady state satisfies the maximal entropy solution, $P(j)=2j$, as it should [1], for any positive $D_{2}$. This FP equation can be written in standard form using the 1st order DC $D_{1}(j)=D(\Delta j)$, as $\partial P/\partial t=(1/2)\partial^{2}D_{2}P/2\partial j^{2}-\partial D_{1}P/\partial j$, where the two DCs are related by $2jD_{1}\left(j\right)=\partial\left(jD_{2}\left(j\right)\right)/\partial j$, as can also be derived by substituting the maximal entropy solution.

Note that the EOM do not contain a true drift term in $j$. However, a “parametric drift” [17] arises because of the gradients in the noise amplitude terms $\boldsymbol{\nu}_{x}(j,\boldsymbol{\Pi})$ (Eqs. 20, 21), which imply that consecutive random steps in the direction of increasing amplitude will tend to be larger than in the opposite direction, even though each individual step is intrinsically symmetric.

3.2 The non-Markovian case

In practice, the Markovian limit is not applicable for the dynamics considered here, since the ACF has a finite correlation time $T_{\mathrm{coh}}$, which can be much longer than the precession period $T_{p}=2\pi/\nu_{p}\left(j\right)$.

The analysis can be extended using several assumptions, which are verified numerically below (Section 4). Since the decay timescale of $C\left(t-t^{\prime}\right)$ is $T_{\mathrm{coh}}$, the response functions mostly contribute to the integral in Eq. (LABEL:eq:P_dot_Nov2) at $|t-t^{\prime}|\ll T_{\mathrm{coh}}<\nu_{j}^{-1}$. Thus, the second term in Eq. (34) is of the order of $\nu_{j}^{2}|t-t^{\prime}|\ll T_{\mathrm{coh}}\nu_{j}^{2}<\nu_{j}$, which is much smaller than the first term ( $\sim\nu_{j}$), and can therefore be neglected. Moreover, on timescales $|t-t^{\prime}|\ll T_{\mathrm{coh}}<\nu_{j}^{-1}$, the state variables $j$, $\phi$ and $u$ (which evolve stochastically on the timescale $\nu_{j}^{-1}$) are almost constant, while $\psi$ evolves deterministically (linearly) in time as $\nu_{p}\left(j\right)(t-t^{\prime})$. Therefore, the response functions can be approximated as

[TABLE]

This gives rise to an effective FP equation

[TABLE]

with an effective DC

[TABLE]

where $S_{\eta}\left(\omega\right)=\mathcal{F}\left[C\left(t\right)\right](\omega)$ is the spectral power density of the noise and $\mathcal{F}$ is the Fourier transform. Since the treatment is only applicable for $t\gg T_{\mathrm{coh}}$, we took the integration limit in Eq. (LABEL:eq:P_dot_Nov2) to infinity.

Note that in the limit where there is no precession i.e. $\nu_{p}\to 0$, the diffusion coefficient approaches the Markovian limit $D_{2}\left(j\right)=2\nu_{j}^{2}\left(j\right)T_{\mathrm{coh}}$. Previous studies (e.g. [10]) identified the diffusion time associated with resonant relaxation as $T_{RR}\sim J_{c}^{2}/\left(\tau_{N}^{2}T_{\mathrm{coh}}\right)\sim 1/D_{2}$, which is equivalent to assuming the Markovian limit. This is a valid approximation in the Newtonian case when $\nu_{GR}/j^{2}<2\pi/T_{\mathrm{coh}}$ 444It is assumed here that mass precession does not induce adiabatic invariance. This remains to be studied in detail..

This effective FP equation has the same form and steady state solution as the one obtained in the Markovian limit, but with a DC that depends on the noise model through its correlation function. Generally, the properties of $D_{2}$ will depend on the “smoothness” of the noise, as reflected by the behavior of the ACF at $t\to 0$ (Figure 1). For example, when the noise is generated by an Ornstein–Uhlenbeck process [24], the noise is continuous but not differentiable, and therefore the ACF is not differentiable at $t=0$, and has the form $C\left(t\right)=e^{-|t|/T_{\mathrm{coh}}}$. The DC is then

[TABLE]

However, such a non-differentiable ACF cannot provide an exact description of the dynamics. The physical noise is generated by the orbital motion of the background stars, and is therefore inherently smooth, with an infinitely differentiable ACF at $t=0$. The smoothness implies that all the derivatives, $\mathrm{d^{n}}C\left(s\right)/\mathrm{d}t^{n}|_{t=0}=\left(2\pi\right)^{-1}\int_{-\infty}^{+\infty}(i\nu)^{n}S_{\eta}(\nu)\mathrm{d}\nu$, are finite (all odd orders are zero by symmetry), so that $S_{\eta}(\nu)$ must fall faster than polynomial. That is, the integrated power of the noise is negligible beyond some characteristic frequency $\nu_{\max}$. When $\nu_{p}>\nu_{\max}$, there is little power in the noise at the precession frequency, and the evolution of $j$ decouples from the background, as is reflected by the diffusion coefficient, $D_{2}\propto S_{\eta}(\nu_{p})=\mathcal{F}\left[C\left(t\right)\right](\nu_{p})$ (Eq. 41). Conversely, the Fourier transform of a non-differentiable noise (with a non-differentiable ACF at $t=0$) contains all frequencies, and in particular has a non-vanishing power at all precession frequencies, which implies a net, noise-driven evolution of $j$ across all phase-space. The sharp suppression of diffusion when $\nu_{p}>\nu_{\max}$ is a manifestation of the general concept of adiabatic invariance, where the action conjugate to a fast dynamical angle is approximately conserved under slow parametric changes of the Hamiltonian [25].

Generally, a physically valid ACF must have the perturbative form of $C\left(t\right)=1-\frac{1}{2}\left\langle\dot{\eta}_{i}^{2}\left(t\right)\right\rangle t^{2}+{\cal O}(t^{4})$, where $\left\langle\dot{\eta}_{i}^{2}\left(t\right)\right\rangle\propto\nu_{\max}^{2}$ is the variability frequency of the noise555For example, for a white noise model truncated at $\nu_{\max}$, $T_{\mathrm{coh}}=\pi/2\nu_{\max}$, $\left\langle\dot{\eta}_{i}^{2}\left(t\right)\right\rangle=\ddot{-C(0)}=\nu_{\max}^{2}/3$., which is not necessarily related to $T_{\mathrm{coh}}$. The angular momentum scale $j_{0}$, where adiabatic invariance becomes important, is therefore given by $\nu_{p}(j_{0})\sim\nu_{\max}$. In the context of relativistic dynamics, and noise models with a single scale $T_{\mathrm{coh}}\propto 1/\nu_{\max}$, which are considered here for simplicity, it is given by $\nu_{GR}/j_{0}^{2}=2\pi/T_{\mathrm{coh}}$, or

[TABLE]

For example, the Gaussian ACF, $C\left(t\right)=\exp\left(-\pi\left(t/T_{\mathrm{coh}}\right)^{2}/4\right)$, where $T_{\mathrm{coh}}=\sqrt{\pi/2\left\langle\dot{\eta}_{i}^{2}\left(t\right)\right\rangle}$, corresponds to a smooth noise with a DC

[TABLE]

which decays very strongly as the precession period falls below the coherence time.

Different noise models and their corresponding ACFs and DCs of the $l=1$ relativistic Hamiltonian $H_{\eta}^{GR}$ (Eq. 17), where $\nu_{p}(j)=\nu_{GR}/j^{2}$ and $\nu_{j}\propto\sqrt{1-j}$, are shown in Figure 1. The results illustrate the relations between the smoothness of the noise, the continuity of the ACF at $t=0$ and the suppression of the DC. The implications of this suppression to the dynamics of relativistic orbits described by $H_{\eta}^{GR}$ are discussed in detail in Section 4.

4 A relativistic induced barrier in phase space

The dynamics of a test star in a stellar cusp around a MBH depend on the dynamics of the background stars, described here in terms of the ACF of a stochastic noise. We analyze here the implications of two specific noise models on the dynamics of the $l=1$ relativistic Hamiltonian (Eq. 17). These represent two limits, that of a continuously differentiable noise, with a Gaussian ACF and that of a non-differentiable noise, with an exponential ACF.

4.1 Time evolution of the cumulative distribution function

Figure 2 shows the evolution of the cumulative distribution function (CDF) at logarithmic time intervals for the two noise models. Time is measured in terms of the global diffusion timescale, $T_{D}=\nu_{j}^{-2}\left(j=0\right)/T_{\mathrm{coh}}$, that is, the inverse of the DC in the Markovian limit. The CDF was obtained by two methods:

Numerical integration of the effective FP equation (40), starting from a narrow PDF (reflecting the initial scatter on timescale $T_{\mathrm{coh}}$, see below) centered around $j_{i}=0.9$ and zero-flux boundary conditions at $j=1$ and at $j_{\min}\to 0$. 2. 2.

Numerical integration of the EOM666Using the vode_f90 [26] implementation of the vode solver [27]. (Eqs. 18, 19) are carried out for many realizations of the noise in a Monte-Carlo fashion. Continuous noise realizations were generated from the ACF using a discrete Fourier transform with randomly drawn phases (e.g. [28]). The simulations were started with an initial $j_{i}=0.9$ and a random orientation of the orbit. Several values of $T_{\mathrm{coh}}\in\left[0.01,1\right]$ and $\nu_{GR}\in\left[0.1,1\right]$ were used to confirm that the evolution of the CDF depends only on $j_{0}=\sqrt{T_{\mathrm{coh}}\nu_{GR}/2\pi}$.

Figure 2 shows that the dynamics of the system crucially depend on the noise model. Under the Gaussian ACF (differentiable) noise, the rate of evolution of the PDF at $j<j_{0}$ is exponentially suppressed since the local diffusion time diverges as $D_{2}^{-1}(j)\propto\exp\left[-4\pi\left(j_{0}/j\right)^{4}\right]$. When the stars are initially placed above $j_{0}$, the PDF reaches a quasi-steady state on a timescale of $T_{D}$, where it drops rapidly to zero at $j_{b}(t)\lesssim j_{0}$. Figure 3 shows that the effective boundary of the distribution $j_{b}\left(t\right)$, defined as the maximum of $dP\left(j,t\right)/dj$, is effectively constant in time, since as we show in A, $j_{b}(t)\approx\left[1+\log\left(t/T_{0}\right)/16\pi\right]^{-1}$ where $j_{b}\left(T_{0}\right)=j_{0}$ and $T_{0}\approx 120j_{0}^{2}T_{D}$. The locus $j_{0}(a)$ can then be interpreted as an effective barrier in phase space, that is, stellar trajectories that cross $j_{0}$ from above, spend only a vanishingly small time below it. In marked contrast, under the exponential ACF (non-differentiable) noise, the PDF evolves much faster and is suppressed only as $D_{2}^{-1}(j)\propto\left(j_{0}/j\right)^{4}$. It does not display a definite barrier, and the fall of the PDF at low-$j$ merely reflects the finite time of the simulation (Figure 2). Note for comparison that the rate of evolution under uncorrelated (discontinuous) noise, equivalent to the case of no GR precession, is nearly uniform in $j$ for $j<1$ (Figure 1), and therefore evolves even faster.

4.2 Evolution of low-$j$ orbits

We integrated the system starting with an initial $j_{i}$ close to $j_{0}$. Figure 4 shows that the evolution of the PDF strongly depends on the noise model and on the location of $j_{i}$. For the Gaussian ACF (differentiable) noise, test stars with initial $j_{i}\gg j_{0}$ reach the equilibrium maximal entropy distribution by time $t\approx T_{D}$, while stars starting at $j_{i}\ll j_{0}$ remain near $j_{i}$, out of equilibrium, on times $t\gg T_{D}$. We model here the evolution from these two opposite regimes and compare results from the FP equation and direct integration of the stochastic EOM. This comparison must however take into account the fact that the FP is applicable only for $t>T_{\mathrm{coh}}$, while the directly integrated orbits evolve on short timescales. We match the initial conditions of the two methods by choosing the initial PDF of the FP equation to have a specific shape and width that reflects the evolution of the integrated orbits up to $t\sim T_{\mathrm{coh}}$. An ensemble of stars with the same initial $j_{i}$ will rapidly spread in $j$ after a short time $t<T_{\mathrm{coh}}$, while the noise is almost constant and $j$ evolves linearly. The evolution of each realization will depend on the specific initial values of the test star’s orbital orientation and the value of the noise . To estimate the PDF at $t=T_{\mathrm{coh}}$, we define $g=j_{i}/j-1$ and use the fact that on timescale $t\ll T_{\mathrm{coh}}$, the Hamiltonian $H_{\eta}^{GR}$ (17) is time independent. Thus

[TABLE]

which implies that $g$ is a random Gaussian with zero mean and $\sigma_{g}^{2}\approx 2j_{i}^{2}\left(\nu_{j}\left(j_{i}\right)/\nu_{GR}\right)^{2}\left(1-\left\langle\cos\left(\psi_{i}-\psi\right)\right\rangle\right)$ where we assume that $j\approx j_{i}$ and $\boldsymbol{\Pi}\approx\boldsymbol{\Pi}_{i}$. For $t\gg j^{2}/\nu_{GR}$, the phase $\psi$ is randomized and $\left\langle\cos\left(\psi_{i}-\psi\right)\right\rangle\approx 0$. For $t\ll j^{2}/\nu_{GR}$, $\psi\approx\psi_{i}+\nu_{GR}t/j^{2}$ and $\left\langle\cos\left(\psi_{i}-\psi\right)\right\rangle\approx 1-\nu_{GR}^{2}t^{2}/\left(2j^{4}\right)$. Therefore we can estimate the $j$ distribution after the coherent phase $t\sim T_{\mathrm{coh}}$ by

[TABLE]

with

[TABLE]

$P_{0}(j;j_{i},\sigma_{g})$ therefore provides the effective initial conditions that correctly match those of the stochastic EOM, as verified by the results shown in Figures 2 and 4.

5 Discussion

The description of the dynamics near a MBH involves many processes, such as two-body relaxation, resonant relaxation, GR-precession and mass precession. Direct relativistic $N$-body simulations include, in principle, all these effects, but are very hard to interpret since these complex dynamics are entangled, and the computational costs limit the simulations to small $N$.

In this study we addressed key elements of this problem analytically, and focused on the interplay between the deterministic GR precession and the stochastic fluctuations of the background torques. We neglected stochastic two-body energy and angular momentum relaxation, since these operate on much longer timescales. We also omitted precession due to the extended stellar mass around the MBH, since it is negligible in the $j\to 0$ limit that is particularly relevant for relativistic orbits.

We demonstrated how the complicated dynamics of a nearly-Keplerian $N$-body system can be described and studied in a formal statistical mechanics framework. We described the dynamics of a test star by a stochastic orbit-averaged Hamiltonian, where the noise terms represent the time-dependent evolution of the background. We extended the usual Legendre multipole expansion in real space to include the new degree of freedom that is introduced by orbit-averaging, namely the orientation of the Keplerian ellipse in its plane, by expanding the Hamiltonian in Wigner matrices. We then derived the explicit stochastic 3D EOM from the first order ($l=1$) relativistic Hamiltonian $H_{\eta}^{GR}$, and presented the corresponding effective FP equation for a general correlated Gaussian noise. We validated that the numerical integration of this FP equation reproduces the statistical properties of the stochastic EOM, and importantly, converges to the maximal entropy limit for $t\to\infty$ (Section 4).

We showed that evolution toward low-$j$ orbits under the $H_{\eta}^{GR}$ Hamiltonian with Gaussian ACF noise (i.e. smooth noise) is restricted by adiabatic invariance. Stars that are initially on high-$j$ orbits where the GR precession time is long, evolve rapidly to lower $j$ and to faster precession, until they reach a threshold at $j_{0}$, where the precession time falls below the coherence time of the background. Beyond that point, the probability of finding a star is vanishingly small, because the fast precession effectively decouples the orbit from the effects of the noise through the mechanism of adiabatic invariance. This is formally expressed by the strong suppression of the diffusion coefficients below $j_{0}$ (steeper than exponential in $j$) (Section 4.1). The phase-space locus $j=j_{0}(a)$ is not a reflecting boundary, but a barrier in the sense the steps toward $j_{0}$ become infinitely small, while steps away from $j_{0}$ become larger. Due to this asymmetry, stars will spend on average only a short time near the barrier. However, we demonstrate that if stars start initially at $j<j_{0}$, for example by being tidally captured there in a binary disruption event, the diffusion timescale back to $j>j_{0}$ is so long that they stall near their capture orbit (Section 4.2). Note that this result does not take into account the fact that stochastic 2-body relaxation, not included in our treatment here, will eventually push the star to higher $j$ and therefore to faster evolution.

This should be contrasted with the behavior of the system under exponential ACF noise (i.e. non-differentiable noise), where there is no barrier, and the deviation from steady state reflects only the finite age of the system (Section 4; Figure 2). Since the system evolves toward isothermal steady state ($P(j)=2j$), where most of the stars are at high-$j$, the trajectories of stars that start out at low-$j$ tend to end up at higher-$j$, while stars that start out at high-$j$, tend to remain there. Such a behavior was observed in the Monte Carlo simulation [7, 29] that used a continuous, but not continuously-differentiable noise model (see [11] for details), which most closely corresponds to our exponential ACF noise model.

The main limitations of this study are that we considered only the $l=1$ term in the Hamiltonian and introduced the physical properties of the noise (smoothness and timescales, i.e. the form of the ACF) as free parameters. However, the maximal variability frequency of the noise, and the phase-space locus of the barrier, can be estimated by general considerations. The noise is a function of the orbital elements of the background stars (Eq. 4), which evolve at a rate that is a combination of the deterministic precession due to GR, and due to the enclosed mass, and a stochastic precession due to the residual torques themselves. A physical correlation function (i.e. differentiable at the origin) can be expanded around $t^{\prime}=0$ by

[TABLE]

Recall that the barrier phenomenon is related to the maximal variability frequency of the noise (Section 3.2). Assume that the fastest precession rate of a typical (i.e. not particularly eccentric) background star is due to mass precession, $\nu_{M}(a)\propto N(a)$ (ignoring GR precession of the background stars). This deterministic in-plane precession is typically faster by $\nu_{M}/\nu_{j}\propto\sqrt{N}$ (Eqs. 17, 22–25) than the stochastic evolution of the other orbital elements [10], and it therefore dominates the noise evolution rate $\sqrt{\left\langle\dot{\eta}^{2}\right\rangle}$. To estimate its magnitude, we note that the torque on a test star with semi-major axis $a$ is dominated by contributions from background stars at $\sim 2a$ [19], and the relevant background precession rate is $\sqrt{\left\langle\dot{\eta}^{2}\right\rangle}\approx\nu_{M}\left(2a,e_{m}\right)$ where here the median eccentricity $e_{m}=\sqrt{1/2}$ was taken as a characteristic value. Under these assumptions, for a Gaussian ACF $T_{\mathrm{coh}}=\sqrt{\pi/2\left\langle\dot{\eta}^{2}\right\rangle}\approx\sqrt{\pi/2}\nu_{M}^{-1}$, and the barrier is located at $j_{b}^{2}=\nu_{GR}\nu_{M}^{-1}/\sqrt{8\pi}\propto 1/\left(N\left(a\right)a\right)$. Note that this scaling is different from a previous empirical determination of $j_{b}\approx 0.6\left(a/1\mathrm{mpc}\right)^{-3/2}$, which was based on a qualitative fit to $N$-body results for an isothermal ($N(a)\propto a$) stellar cusp model [11, 12]. For that system, our analysis777For an $\alpha=2$ cusp with noise described by a Gaussian ACF and assuming $T_{\mathrm{coh}}=\sqrt{\pi/2}\nu_{M}^{-1}\left(2a,\sqrt{1/2}\right)$ , $j_{b}^{2}\simeq 0.4(r_{g}/a)(Q/N(a))$. Alternatively, assuming $T_{\mathrm{coh}}=Q\nu_{r}^{-1}\left(a\right)/\sqrt{N(2a)}$, $j_{b}^{2}\simeq 0.33\left(r_{g}/a\right)\left(Q/\sqrt{N\left(a\right)}\right)$ indicates that the barrier should be at $j_{b}\approx 0.6\left(a/1\mathrm{mpc}\right)^{-1}$. Alternatively, if the coherence time of the noise reflects the stochastic evolution of the background stars due to the residual torque themselves, then $T_{\mathrm{coh}}\approx Q\nu_{r}^{-1}\left(a\right)/\sqrt{N(2a)}$, and an even flatter relation $j_{b}\approx 0.4a^{-3/4}$ is predicted. Either loci appear to be reasonably consistent with the presently available numeric data, which has a limited range in $a$ and $j$. Recently, Hamers et al. [30], measured the diffusion coefficient from a suite of simulations. Based on these diffusion coefficients, they presented two expression for the barrier location. One is consistent with [11] and the other is consistent with the analysis presented here. However, using the available data they could not rule out either one. Therefore, formal fits and better statistics will be required to discriminate between the different models.

It should be emphasized that the neglected higher ($l>1$) terms in the multipole expansion of the Hamiltonian, while smaller in magnitude, could in principle introduce shorter timescales in the noise, and correspondingly result in diffusion coefficients that decay at smaller $j$. The superposition of many such decaying terms could result in a slower, power-law decay, which will blur the barrier. Further study is required to apply the insights obtained here about dynamics driven by correlated noise, and in particular the close relation between the smoothness of the background noise, adiabatic invariance, and very low-angular momentum orbits, to real physical systems.

We acknowledge support by the ERC Starting Grant No. 202996, DIP-BMBF Grant No. 71-0460-0101, under which this project was initiated, and the I-CORE Program of the PBC and ISF (Center No. 1829/12).

Appendix A Time evolution of the barrier for a Gaussian ACF noise

The location of the barrier $j_{b}\left(t\right)$ can be defined as the maximum of $\partial P\left(j,t\right)/\partial j$. Assuming that for $j\approx j_{b}$, $P\left(j,t\right)$ is a self-similar function of $j-j_{b}$, that is $P\left(j,t\right)\approx P\left(j-j_{b}\left(t\right)\right)$, and $\int_{j_{b}}^{1}P\left(j,t\right)dj$ is constant in time, we obtain

[TABLE]

and

[TABLE]

Using Eq. (40) and $\partial^{2}P\left(j_{b}\left(t\right),t\right)/\partial j^{2}=0$ we obtain

[TABLE]

For the Gaussian ACF assuming $\nu_{p}=\nu_{GR}/j^{2}$ the DC is given by (Eq. 44),

[TABLE]

where $T_{D}=\nu_{j}^{-2}\left(j=0\right)/T_{\mathrm{coh}}$. By defining $x=j_{0}/j_{b}$ and $s=t/T_{D}$, we use Eq. (51) to obtain

[TABLE]

where we assumed $x\approx 1$. Therefore, as demonsrated in Figure 3, the evolution of the barrier can be approximated by a logarithmically-suppressed evolution function

[TABLE]

where

[TABLE]

is the time the barrier reaches the point $j_{b}\left(t\right)=j_{0}$. Note that since $j_{0}^{2}T_{D}=\nu_{GR}\nu_{j}^{-2}\left(j=0\right)/2\pi$, $T_{0}$ is independent of $T_{\mathrm{coh}}$.