Vector Resonant Relaxation of Stars around a Massive Black Hole

Most nearby galaxies possess a massive black hole (MBH) in their center, surrounded by a nuclear stellar cluster (NSC; Genzel et al. 2010; Kormendy & Ho 2013; Graham 2016). The dynamical evolution of the stellar cluster comprises numerous processes acting on different timescales (Rauch & Tremaine 1996; Hopman & Alexander 2006; Alexander 2017; see also Figure 1 in Kocsis & Tremaine 2011). Since the gravitational potential is dominated by the central MBH, stars move on nearly Keplerian orbits. The deviations from a Keplerian potential due to the stellar potential and the relativistic corrections cause the Keplerian ellipses to precess in their orbital plane. Subsequently, through the non-spherical components of the potential fluctuations, the orbital orientations of the stars get reshuffled, without changing the magnitude of their angular momentum or their Keplerian energy, through a process called vector resonant relaxation (VRR; Kocsis & Tremaine 2015 and references therein), which is the focus of this work. Resonant torques' coupling between the precessing stars then leads to a diffusion of the stars' angular momentum magnitude, a process called scalar resonant relaxation (SRR; Bar-Or & Fouvry 2018 and references therein). Finally, on longer timescales, close encounters between stars lead to the relaxation of their Keplerian energy and angular momentum (Bahcall & Wolf 1976, 1977; Lightman & Shapiro 1977; Cohn & Kulsrud 1978; Shapiro & Marchant 1978).

VRR can be a driving force behind several dynamical phenomena in galactic centers, including the warping of accretion (Bregman & Alexander 2009, 2012) and stellar (Kocsis & Tremaine 2011) disks, as well as a catalyzer of binary mergers (Hamers et al. 2018). Since its first presentation by Rauch & Tremaine (1996), VRR was studied numerically using both full (Eilon et al. 2009) and orbit-averaged (Kocsis & Tremaine 2015) numerical simulations. More recently, Roupas et al. (2017), Takács & Kocsis (2018), and Szölgyén & Kocsis (2018) studied the thermodynamical equilibria of VRR, and Fouvry et al. (2019) its axisymmetric limit.

In this paper, building upon these works, we set out to offer a detailed characterization of the VRR process in the limit of an isotropic distribution of stars. To do so, in Section 2, we present the fundamental equations of VRR. In Section 3, we characterize the properties of the potential fluctuations in the system, as inferred from estimates of the correlation function at the initial time. This will allow us to describe in Section 4 the random walk of a test particle's orientation, and develop an effective stochastic equation of motion which can efficiently mimic these random motions. Detailed comparisons of these results with numerical simulations are presented throughout these sections. In Section 5, we detail the important self-consistency existing between the potential fluctuations and the properties of the orientations' random walks. In Section 6, we use this new formalism to present the timescales associated with VRR in an NSC similar to the Milky Way's. We conclude in Section 7.

We consider a set of N stars orbiting an MBH of mass ${M}_{\bullet }$. On timescales longer than the in-plane precession but shorter than the time to change the orbital eccentricity (by SRR) and the time to change the semimajor axis (by two-body relaxation) the mutual interactions between two stars can be orbit-averaged over their respective mean anomalies and in-plane precession angles. As a result, each star can be replaced by a disk of mass m extending between ${r}_{{\rm{p}}}$ = (1−e)a and ${r}_{{\rm{a}}}$ = (1 + e)a with surface density ${\rm{\Sigma }}(r)$ = ${[2{\pi }^{2}a\sqrt{r-{r}_{{\rm{p}}}}\sqrt{{r}_{{\rm{a}}}-r}]}^{-1}$, where the semimajor axis a and eccentricity e can be assumed to be constant in time (see the illustration in Figure 1). Following this double orbit-average, one can associate to each star a set of conserved quantities ${\boldsymbol{K}}$ = (m, a, e) and a time-dependent normal vector $\widehat{{\boldsymbol{L}}}$, with ${\boldsymbol{L}}=L\widehat{{\boldsymbol{L}}}$ the orbital angular momentum. We introduce the spherical coordinates (θ, ϕ), so that $\widehat{{\boldsymbol{L}}}\,=(\sqrt{1-{u}^{2}}\cos \phi ,\sqrt{1-{u}^{2}}\sin \phi ,u)$, with $u=\cos (\theta )$. Studying VRR then amounts to studying the long-term dynamics of each star's normal vector $\widehat{{\boldsymbol{L}}}$.

Zoom In Zoom Out Reset image size

Following Kocsis & Tremaine (2015), the effective single-particle Hamiltonian of VRR reads

$$\begin{eqnarray}\begin{array}{rcl}H & = & \displaystyle \sum _{i=1}^{N}{\left\langle -\displaystyle \frac{{{Gmm}}_{i}}{| {\boldsymbol{r}}(t)-{{\boldsymbol{r}}}_{i}({t}^{{\prime} })| }\right\rangle }_{t,{t}^{{\prime} }}\\ & = & \ -L({\boldsymbol{K}})\displaystyle \sum _{\mathop{{\ell }\geqslant 2}\limits_{\mathrm{even}}}\displaystyle \sum _{m=-{\ell }}^{{\ell }}{M}_{{\ell }m}({\boldsymbol{K}},t)\,{Y}_{{\ell }m}(\widehat{{\boldsymbol{L}}})\end{array}\end{eqnarray} \tag{ 1 }$$

with ${\boldsymbol{r}}$(t), ${\boldsymbol{r}}$_i(t) the positions of the test star and the star i as they move along their (in-plane) precessing Keplerian orbits, $\langle \,\cdot \,{\rangle }_{t,{t}^{{\prime} }}$ the double orbit-average over these motions, and ${\boldsymbol{K}}$ and K_i their respective conserved parameters. In the second line of Equation (1), we introduced the magnetizations

$$\begin{eqnarray}&&{M}_{{\ell }m}({\boldsymbol{K}},t)=\displaystyle \sum _{i=1}^{N}{{ \mathcal J }}_{{\ell }}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}_{i}\right]\,{Y}_{{\ell }m}({\widehat{{\boldsymbol{L}}}}_{i}(t)),\end{eqnarray} \tag{ 2 }$$

where the coupling coefficients ${{ \mathcal J }}_{{\ell }}$[${\boldsymbol{K}}$, K_i] are defined in Equation (51) below, and we used real spherical harmonics Y_ℓm($\widehat{{\boldsymbol{L}}}$) (defined in Equation (53)). It is also important to note that only even harmonics with ℓ ≥ 2 contribute to the particles' dynamics, by virtue of the symmetries of the interaction.

Hamilton's equations of motion read

$$\begin{eqnarray}&&\displaystyle \frac{d\phi }{dt}=\displaystyle \frac{{\rm{\partial }}H}{{\rm{\partial }}{L}_{z}};\,\displaystyle \frac{d{L}_{z}}{dt}=-\displaystyle \frac{{\rm{\partial }}H}{{\rm{\partial }}\phi },\end{eqnarray} \tag{ 3 }$$

where L_z = Lu is an action and ϕ its conjugated angle. The evolution of the angular momentum of a single test particle is given by

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{L}}}{{dt}}=\displaystyle \frac{\partial }{\partial {\boldsymbol{L}}}\times (H{\boldsymbol{L}})=L\displaystyle \sum _{{\ell },m}{M}_{{\ell }m}({\boldsymbol{K}},t)\,{{\boldsymbol{X}}}_{{\ell }m}(\widehat{{\boldsymbol{L}}}(t)),\end{eqnarray} \tag{ 4 }$$

where the first equality is a compact rewriting of the individual equations of motion. In the second equality, we also introduced ${\boldsymbol{X}}$_ℓm($\widehat{{\boldsymbol{L}}}$) = $\widehat{{\boldsymbol{L}}}$ × ∂Y_ℓm($\widehat{{\boldsymbol{L}}}$)/∂$\widehat{{\boldsymbol{L}}}$ as the real vector spherical harmonics. Because L is constant, one has $d\widehat{{\boldsymbol{L}}}/{dt}={L}^{-1}d{\boldsymbol{L}}/{dt}$.

Inspired by Klimontovich (1967), the state of the system of N stars at time t is fully characterized by the discrete distribution function (DF) 

$$\begin{eqnarray}&&\varphi (\widehat{{\boldsymbol{L}}},{\boldsymbol{K}},t)=\displaystyle \sum _{i}{\delta }_{{\rm{D}}}(\widehat{{\boldsymbol{L}}}-{\widehat{{\boldsymbol{L}}}}_{i}(t))\,{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}_{i}),\end{eqnarray} \tag{ 5 }$$

with ${\delta }_{{\rm{D}}}$ the Dirac delta, ${\delta }_{{\rm{D}}}(\widehat{{\boldsymbol{L}}}-{\widehat{{\boldsymbol{L}}}}_{i})$ = ${\delta }_{{\rm{D}}}(u-{u}_{i}){\delta }_{{\rm{D}}}(\phi -{\phi }_{i})$, and ${\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}_{i})$ = ${\delta }_{{\rm{D}}}(m-{m}_{i}){\delta }_{{\rm{D}}}(a-{a}_{i}){\delta }_{{\rm{D}}}(e-{e}_{i})$, with the associated volumes $d\widehat{{\boldsymbol{L}}}$ = $d$u $d$ϕ and $d$ ${\boldsymbol{K}}$ = $d$m $d$a $d$e. The continuity equation, ${\rm{\partial }}\varphi /{\rm{\partial }}t$ = $-{\rm{\partial }}/{\rm{\partial }}\hat{{\boldsymbol{L}}}\cdot [\varphi \,d\hat{{\boldsymbol{L}}}/dt]$, then gives us

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial \varphi (\widehat{{\boldsymbol{L}}},{\boldsymbol{K}},t)}{\partial t} & = & -\displaystyle \sum _{{\ell },m}\displaystyle \int d{\widehat{{\boldsymbol{L}}}}^{{\prime} }{d}{{\boldsymbol{K}}}^{{\prime} }\,\varphi ({\widehat{{\boldsymbol{L}}}}^{{\prime} },{{\boldsymbol{K}}}^{{\prime} },t)\,{{ \mathcal J }}_{{\ell }}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} }\right]\\ & & \times \ {Y}_{{\ell }m}({\widehat{{\boldsymbol{L}}}}^{{\prime} })\,{{\boldsymbol{X}}}_{{\ell }m}(\widehat{{\boldsymbol{L}}})\cdot \displaystyle \frac{\partial \varphi (\widehat{{\boldsymbol{L}}},{\boldsymbol{K}},t)}{\partial \widehat{{\boldsymbol{L}}}},\end{array}\end{eqnarray} \tag{ 6 }$$

where we use the fact that the vector spherical harmonics satisfy ∂/∂$\widehat{{\boldsymbol{L}}}$ · ${\boldsymbol{X}}$($\widehat{{\boldsymbol{L}}}$) = 0. This equation can subsequently be developed in spherical harmonics by writing

$$\begin{eqnarray}&&\varphi (\widehat{{\boldsymbol{L}}},{\boldsymbol{K}},t)={\varphi }_{\alpha }({\boldsymbol{K}},t)\,{Y}_{\alpha }(\widehat{{\boldsymbol{L}}}),\end{eqnarray} \tag{ 7 }$$

where the sum over the index α = (ℓ_α, m_α) is implied, and we introduce

$$\begin{eqnarray}&&{\varphi }_{\alpha }({\boldsymbol{K}},t)=\displaystyle \sum _{i}{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}_{i})\,{Y}_{\alpha }({\widehat{{\boldsymbol{L}}}}_{i}(t)).\end{eqnarray} \tag{ 8 }$$

When expanded in spherical harmonics, Equation (6) becomes

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\varphi }_{\alpha }({\boldsymbol{K}},t)}{\partial t}=-\int d{{\boldsymbol{K}}}^{{\prime} }\,{Q}_{\alpha \gamma \delta }({\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} })\,{\varphi }_{\gamma }({{\boldsymbol{K}}}^{{\prime} },t)\,{\varphi }_{\delta }({\boldsymbol{K}},t),\end{eqnarray} \tag{ 9 }$$

where the sums over the harmonic indices (γ, δ) are implied, and we introduce the time-independent coupling tensor Q_αγδ(${\boldsymbol{K}}$, ${{\boldsymbol{K}}}^{{\prime} }$) = ${{ \mathcal J }}_{{{\ell }}_{\gamma }}$[${\boldsymbol{K}}$, ${{\boldsymbol{K}}}^{{\prime} }$] E_αγδ, with E_αγδ the (real) Elsässer coefficients (James 1973) (see Appendix B for their properties)

$$\begin{eqnarray}&&{E}_{\alpha \gamma \delta }=\int d\widehat{{\boldsymbol{L}}}\,{Y}_{\alpha }(\widehat{{\boldsymbol{L}}})\,{{\boldsymbol{X}}}_{\gamma }(\widehat{{\boldsymbol{L}}})\cdot \displaystyle \frac{\partial {Y}_{\delta }(\widehat{{\boldsymbol{L}}})}{\partial \widehat{{\boldsymbol{L}}}}.\end{eqnarray} \tag{ 10 }$$

Equation (9) is an exact writing of the fundamental evolution equation for VRR. Its complexity stems in particular from being a quadratic matrix differential equation in the fields ${\varphi }_{\alpha }({\boldsymbol{K}},t)$.

All the upcoming derivations will be illustrated by comparisons with numerical simulations. In Appendix C, we present the fiducial system considered, as well as the details of our numerical implementation. In Figure 2, we illustrate a subset of trajectories from one such simulation.

Zoom In Zoom Out Reset image size

As a first step toward the characterization of the correlated stochastic dynamics of one star in that system, we focus on describing the properties of the density fluctuations generated as a whole by the system's N particles. In particular, we will show how one can use estimates of the derivatives of the correlation function of the density fluctuations at the initial time to provide a sensible ansatz (see Equation (22)) for the time dependences of this same correlation function.

The harmonic coefficients ${\varphi }_{\alpha }({\boldsymbol{K}},t)$ (Equation (8)) describe the full state of the N-particle system (N ≫ 1) at time t and can therefore be treated as stochastic density fluctuations, assumed to be Gaussian random fields. Assuming that the system's evolution is stationary in time, the properties of these fluctuations are captured by the correlation function

$$\begin{eqnarray}&&{C}_{\alpha \beta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,t-{t}^{{\prime} })\equiv \left\langle {\varphi }_{\alpha }({\boldsymbol{K}},t)\,{\varphi }_{\beta }({\boldsymbol{K}}^{\prime} ,{t}^{{\prime} })\right\rangle ,\end{eqnarray} \tag{ 11 }$$

where $\left\langle \cdot \right\rangle$ is the ensemble average over realizations (initial conditions and trajectories of the N particles).

The correlation function is even with respect to time and generically decreases to zero on a timescale larger than some coherence time ${T}_{{\rm{c}}}$. As a result, as a first approximation, it is therefore reasonable to replace C_αβ(${\boldsymbol{K}}$, ${\boldsymbol{K}}$', t − ${t}^{{\prime} }$) by a Gaussian function, tailored to match the function's behavior for t ≪ ${T}_{{\rm{c}}}$; see Figure 3 for a justification.

Zoom In Zoom Out Reset image size

In Appendix D, we compute the first two derivatives of the correlation function, and we show in Equations (64) and (68) that

$$\begin{eqnarray}&&{C}_{\alpha \beta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,0)={\delta }_{\alpha }^{\beta }\,{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{\boldsymbol{K}}^{\prime} )\,n({\boldsymbol{K}})\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}{C}_{\alpha \beta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,t)\right|}_{t=0}=-{\delta }_{\alpha }^{\beta }\,{A}_{{{\ell }}_{\alpha }}{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{\boldsymbol{K}}^{\prime} )\,n({\boldsymbol{K}})\,{{\rm{\Gamma }}}^{2}({\boldsymbol{K}}),\end{eqnarray} \tag{ 13 }$$

with the coefficient

$$\begin{eqnarray}&&{A}_{{\ell }}={\ell }({\ell }+1).\end{eqnarray} \tag{ 14 }$$

In Equation (13), we introduced n(${\boldsymbol{K}}$), the DF of the stars' ${\boldsymbol{K}}$ = (a, e, m) parameters, that satisfies the normalization convention $\int d\widehat{{\boldsymbol{L}}}d{\boldsymbol{K}}n({\boldsymbol{K}})=N$. We also introduced the decay rate of the correlation function Γ(${\boldsymbol{K}}$) as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}^{2}({\boldsymbol{K}})=\int d{\boldsymbol{K}}^{\prime} \,n({{\boldsymbol{K}}}^{{\prime} })\,\displaystyle \sum _{{\ell }}{B}_{{\ell }}\,{{ \mathcal J }}_{{\ell }}^{2}[{\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} }],\end{eqnarray} \tag{ 15 }$$

with the coefficient

$$\begin{eqnarray}&&{B}_{{\ell }}=\displaystyle \frac{{\ell }({\ell }+1)(2{\ell }+1)}{8\pi }.\end{eqnarray} \tag{ 16 }$$

Gathering Equations (12) and (13), we can approximate the correlation function ${C}_{\alpha \beta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,t)$ by

$$\begin{eqnarray}&&{C}_{\alpha \beta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,t)={\delta }_{\alpha }^{\beta }\,{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}^{{\prime} })\,{C}_{{{\ell }}_{\alpha }}({\boldsymbol{K}},t),\end{eqnarray} \tag{ 17 }$$

where C_ℓ(${\boldsymbol{K}}$, t) is a function decaying like a Gaussian:

$$\begin{eqnarray}&&{C}_{{\ell }}({\boldsymbol{K}},t)=n({\boldsymbol{K}})\,{e}^{-\tfrac{{A}_{{\ell }}}{2}{\left(t/{T}_{{\rm{c}}}({\boldsymbol{K}})\right)}^{2}},\end{eqnarray} \tag{ 18 }$$

where we introduce the torque time

$$\begin{eqnarray}&&{T}_{{\rm{c}}}({\boldsymbol{K}})=\displaystyle \frac{1}{{\rm{\Gamma }}({\boldsymbol{K}})}.\end{eqnarray} \tag{ 19 }$$

Equations (17) and (18) are the main result of this section, as they provide us with a simple estimate for the time evolution of the ensemble-averaged correlation function of the fluctuations in the system. In Figure 3 we compare this estimate to the correlations measured in the numerical simulations and, to shorten the main text, we detail the procedure followed to obtain that figure in Appendix H.1. As expected, this estimation matches the numerical measurements on short timescales. Capturing the late-time, non-Gaussian behavior of the correlation function requires a self-consistent determination of the time-dependence of the noise. This is investigated in Section 5, and allows for an improved noise prediction in Figure 3. In Section 4, we will use the previous correlation functions as source terms to describe the dynamics of a test particle embedded in that noisy environment. However, in that section we will see that the ensemble-averaged correlation function does not capture the full dynamics induced on a test particle. Indeed, globally conserved quantities (such as the total energy) prevent the system from being fully ergodic: even after long times the system will not explore the entire realization space and therefore time averages are not equivalent to ensemble averages. For a given realization "${\rm{r}}$," we therefore define the time-averaged correlation function

$$\begin{eqnarray}&&{C}_{\alpha \beta }^{{\rm{r}}}({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,t-{t}^{{\prime} })\equiv {\left\langle {\varphi }_{\alpha }({\boldsymbol{K}},t){\varphi }_{\beta }({\boldsymbol{K}}^{\prime} ,{t}^{{\prime} })\right\rangle }_{T},\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&\langle f{\rangle }_{T}\equiv \displaystyle \frac{1}{T}{\int }_{0}^{T}{dt}\,f(t)\end{eqnarray} \tag{ 21 }$$

stands for the time average over some long timescale T. As previously, we will assume that ${C}_{\alpha \beta }^{{\rm{r}}}$(${\boldsymbol{K}}$, ${\boldsymbol{K}}$', t) can be replaced by

$$\begin{eqnarray}&&{C}_{\alpha \beta }^{{\rm{r}}}({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,t)={\delta }_{\alpha }^{\beta }\,{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{\boldsymbol{K}}^{\prime} )\,{C}_{{{\ell }}_{\alpha }}^{{\rm{r}}}({\boldsymbol{K}},t),\end{eqnarray} \tag{ 22 }$$

with the Gaussian time dependence

$$\begin{eqnarray}&&{C}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}},t)={n}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}})\,{e}^{-\tfrac{{A}_{{\ell }}}{2}{\left(t/{T}_{{\rm{c}}}({\boldsymbol{K}})\right)}^{2}}.\end{eqnarray} \tag{ 23 }$$

Here, we define the effective amplitude ${n}_{{\ell }}^{{\rm{r}}}$(${\boldsymbol{K}}$) as the mean value over m_α, so that

$$\begin{eqnarray}&&{n}_{{{\ell }}_{\alpha }}^{{\rm{r}}}({\boldsymbol{K}})=\displaystyle \frac{1}{2{{\ell }}_{\alpha }+1}\displaystyle \sum _{{m}_{\alpha },\beta }\displaystyle \int d{{\boldsymbol{K}}}_{\beta }\,{C}_{\alpha \beta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },0),\end{eqnarray} \tag{ 24 }$$

and from Equations (12), we have

$$\begin{eqnarray}&&\left\langle {n}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}})\right\rangle =n({\boldsymbol{K}}),\end{eqnarray} \tag{ 25 }$$

which is independent of the considered harmonic. It is important to note that ${n}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}})$ varies between different realizations. In Equation (98), we illustrate how one can compute its variance, and show how this originates from the constraint of total energy conservation.

In Equation (23), we assumed, for simplicity, that the torque time, ${T}_{{\rm{c}}}$(${\boldsymbol{K}}$), is independent of the considered realization. These various choices ensure that the ansatz from Equation (22) satisfies the constraints from Equations (12) and (13) when ensemble-averaged. As highlighted by Equation (22), this correlation is diagonal both w.r.t. the harmonic indices (via ${\delta }_{\alpha }^{\beta }$) and w.r.t. the considered parameters (via ${\delta }_{{\rm{D}}}$(${\boldsymbol{K}}$ − ${\boldsymbol{K}}$')). Following Equation (23), we note that the time-dependence of this correlation is controlled by both an effective amplitude, ${n}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}})$, and a torque time, ${T}_{{\rm{c}}}$(${\boldsymbol{K}}$), which depend on the considered parameter ${\boldsymbol{K}}$.

In the previous section, we characterized the noise fluctuations resulting from the coupled motions of the system's N particles. Assuming that the statistics of this noise follows the correlation function obtained in Equation (17), our goal is now to investigate the stochastic dynamics of one given test particle embedded in that fluctuating environment. In Figure 4, we illustrate one such random walk by highlighting the time evolution of the orientation of a single particle in one fiducial simulation.

Zoom In Zoom Out Reset image size

Throughout this section, we use the test particle limit, i.e., we assume that the motion of the test particle is fully determined by the time-dependent density of the background particles and we neglect any backreaction of the test particle onto the background particles. We denote the parameters of the test particle with ${{\boldsymbol{K}}}_{{\rm{t}}}$, and its orientation at time t with ${\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}$(t). Similarly to Equation (5), the current orientation of the test particle is fully characterized by the single-particle DF: 

$$\begin{eqnarray}&&{\varphi }^{{\rm{t}}}(\widehat{{\boldsymbol{L}}},t)={\delta }_{{\rm{D}}}(\widehat{{\boldsymbol{L}}}-{\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}(t)),\end{eqnarray} \tag{ 26 }$$

which can be expanded as ${\varphi }^{{\rm{t}}}$($\widehat{{\boldsymbol{L}}}$, t) = Y_α($\widehat{{\boldsymbol{L}}}$) ${\varphi }_{\alpha }^{{\rm{t}}}(t)$ (with the sum over α implied), where

$$\begin{eqnarray}&&{\varphi }_{\alpha }^{{\rm{t}}}(t)={Y}_{\alpha }({\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}(t)).\end{eqnarray} \tag{ 27 }$$

In practice, because the first spherical harmonics ℓ = 1 is such that ${\varphi }_{(1,m)}^{{\rm{t}}}$ (t) ∼ ${\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}$(t), characterizing the random walk of ${\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}$, the test particle's orbital orientation, as illustrated in Figure 4, therefore requires knowledge of the correlation properties of ${\varphi }_{\alpha }^{{\rm{t}}}(t)$ for ℓ_α = 1.

The evolution equation for ${\varphi }_{\alpha }^{{\rm{t}}}(t)$ follows from Equation (9), and reads

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {\varphi }_{\alpha }^{{\rm{t}}}(t)}{\partial t} & = & -\displaystyle \int d{\boldsymbol{K}}\,{Q}_{\alpha \gamma \delta }({{\boldsymbol{K}}}_{{\rm{t}}},{\boldsymbol{K}})\,{\varphi }_{\gamma }({\boldsymbol{K}},t)\,{\varphi }_{\delta }^{{\rm{t}}}(t)\\ & = & -{Q}_{\alpha \delta }^{{\rm{t}}}(t)\,{\varphi }_{\delta }^{{\rm{t}}}(t),\end{array}\end{eqnarray} \tag{ 28 }$$

where ${\varphi }_{\gamma }({\boldsymbol{K}},t)$ is the harmonic coefficients of the background particles' field, as defined in Equation (8), and we introduced ${Q}_{\alpha \delta }^{{\rm{t}}}(t)$ = $\int d{\boldsymbol{K}}\,{Q}_{\alpha \gamma \delta }({{\boldsymbol{K}}}_{{\rm{t}}},{\boldsymbol{K}})\,{\varphi }_{\gamma }({\boldsymbol{K}},t)$, which is a time-dependent external forcing term driving the dynamics of the test particle. The time-dependence of this source only originates from the background particles (via ${\varphi }_{\gamma }({\boldsymbol{K}},t)$), whose correlation properties were investigated in the previous section. We also emphasize that this forcing term also depends on ${{\boldsymbol{K}}}_{{\rm{t}}}$, the parameters of the considered test particle.

Equation (28) takes the form of a time-dependent linear matrix differential equation for the test particle's harmonic coefficients ${\varphi }_{\alpha }^{{\rm{t}}}(t)$. In order to guarantee well-behaved asymptotics of the test particle's motion for large times, we approach the resolution of Equation (28) via the Magnus series (see Blanes et al. 2009 for a review). In that framework, one can generically solve for the motion of the test particle as

$$\begin{eqnarray}&&{\varphi }_{\alpha }^{{\rm{t}}}(t)={\left[{e}^{{\rm{\Omega }}({t}^{{\prime} },t)}\right]}_{\alpha \delta }\,{\varphi }_{\delta }^{{\rm{t}}}({t}^{{\prime} }),\end{eqnarray} \tag{ 29 }$$

where the matrix Ω (${t}^{{\prime} }$, t) is constructed as a series expansion of the form ${\rm{\Omega }}({t}^{{\prime} },t)={\sum }_{k\geqslant 1}{{\rm{\Omega }}}_{k}({t}^{{\prime} },t)$, whose first terms are

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Omega }}}_{1}({t}^{{\prime} },t) & = & -{\displaystyle \int }_{{t}^{{\prime} }}^{t}{ds}\,{Q}^{{\rm{t}}}(s),\\ {{\rm{\Omega }}}_{2}({t}^{{\prime} },t) & = & \displaystyle \frac{1}{2}{\displaystyle \int }_{{t}^{{\prime} }}^{t}{{ds}}_{1}{\displaystyle \int }_{{t}^{{\prime} }}^{{s}_{1}}{{ds}}_{2}\left[{Q}^{{\rm{t}}}({s}_{1}),{Q}^{{\rm{t}}}({s}_{2})\right],\end{array}\end{eqnarray} \tag{ 30 }$$

where [A, B] = AB − BA is the matrix commutator. As in Equation (20), for a given realization, the statistics of the motion of the test particle are captured by the stationary time-averaged correlation function

$$\begin{eqnarray}&&{C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}(t-{t}^{{\prime} })\equiv {\left\langle {\varphi }_{\alpha }^{{\rm{t}}}(t){\varphi }_{\beta }^{{\rm{t}}}({t}^{{\prime} })\right\rangle }_{T}.\end{eqnarray} \tag{ 31 }$$

Using Equation (29), we can write this correlation function as

$$\begin{eqnarray}&&{C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}(t)={\left\langle {\left[{e}^{{\rm{\Omega }}(t)}\right]}_{\alpha \delta }\right\rangle }_{T}\,{\left\langle {\varphi }_{\delta }^{{\rm{t}}}(0){\varphi }_{\beta }^{{\rm{t}}}(0)\right\rangle }_{{\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}},\end{eqnarray} \tag{ 32 }$$

where we rely on our test particle's assumption (i.e., independence hypothesis; Corrsin 1959), which allowed us to separate the time average (denoted as $\langle \cdot {\rangle }_{T}$) over the background particles generating the noise, and the average over the initial location of the test particle (denoted as $\langle \cdot {\rangle }_{{\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}}$). We also rely on the hypothesis that the noise is stationary in time, so that $\langle {e}^{{\rm{\Omega }}({t}^{{\prime} },t)}{\rangle }_{T}=\langle {e}^{{\rm{\Omega }}(t-{t}^{{\prime} })}{\rangle }_{T}$, with Ω (t) ≡ Ω (0, t).

In Appendix E, we rely on the cumulant theorem to compute the two averages appearing in Equation (32). This allows us to rewrite the test particle's correlation function as

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}(t) & = & \displaystyle \frac{{\delta }_{\alpha }^{\beta }}{4\pi }\,\exp \left\{-{A}_{{{\ell }}_{\alpha }}\displaystyle \sum _{{\ell }}{B}_{{\ell }}\displaystyle \int d{\boldsymbol{K}}\,{n}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}}){{ \mathcal J }}_{{\ell }}^{2}\left[{{\boldsymbol{K}}}_{{\rm{t}}},{\boldsymbol{K}}\right]\right.\\ & & \times \ \left.\displaystyle \frac{{T}_{{\rm{c}}}^{2}({\boldsymbol{K}})}{{A}_{{\ell }}}\chi \left[\displaystyle \frac{t\sqrt{{A}_{{\ell }}/2}}{{T}_{{\rm{c}}}({\boldsymbol{K}})}\right]\right\}\end{array}\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray*}\begin{array}{r}\begin{array}{c}\simeq \displaystyle \frac{{\delta }_{\alpha }^{\beta }}{4\pi }\times \left\{\begin{array}{cc}{e}^{-{\textstyle {\textstyle \tfrac{{A}_{{\ell }_{\alpha }}}{2}}}\left(\displaystyle \sum _{\ell }{B}_{\ell }\displaystyle \int d{\boldsymbol{K}}{n}_{\ell }^{{\rm{r}}}({\boldsymbol{K}}){J}_{\ell }^{2}[{{\boldsymbol{K}}}_{{\rm{t}}},{\boldsymbol{K}}]\right){t}^{2}} & t\ll {T}_{{\rm{c}}},\\ {e}^{-{\textstyle {\textstyle \tfrac{\sqrt{\pi }{A}_{{\ell }_{\alpha }}}{2}}}\left(\displaystyle \sum _{\ell }\displaystyle \frac{{B}_{\ell }}{\sqrt{{A}_{\ell }/2}}\displaystyle \int d{\boldsymbol{K}}{n}_{\ell }^{{\rm{r}}}({\boldsymbol{K}}){J}_{\ell }^{2}[{{\boldsymbol{K}}}_{{\rm{t}}},{\boldsymbol{K}}]{T}_{{\rm{c}}}({\boldsymbol{K}})\right)t} & t\gg {T}_{{\rm{c}}},\end{array}\right.\end{array}\end{array}\end{eqnarray*}$$

where we introduce the dimensionless function

$$\begin{eqnarray}\begin{array}{rcl}\chi (\tau ) & = & {\displaystyle \int }_{0}^{\tau }d{\tau }_{1}{\displaystyle \int }_{0}^{\tau }d{\tau }_{2}\,{e}^{-{\left({\tau }_{1}-{\tau }_{2}\right)}^{2}}\\ & = & \ -1+{e}^{-{\tau }^{2}}+\tau \,\sqrt{\pi }\,\mathrm{erf}(\,\tau )\\ & \simeq & \ \left\{\begin{array}{cc}{\tau }^{2} & \tau \ll 1,\\ \sqrt{\pi }\tau & \tau \gg 1.\end{array}\right.\end{array}\end{eqnarray} \tag{ 34 }$$

As can be seen from the time-dependence of the exponent in Equation (33), one can note that, on short timescales, t ≪ ${T}_{{\rm{c}}}$, the correlation ${C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}$ decays like a Gaussian and the motion of the particle is ballistic (e.g., (Δ$\widehat{{\boldsymbol{L}}}$)² ∝ t²). On these short timescales, the random walk of the test particle is analogous to that induced by a time-independent fluctuation. On long timescales, t ≫ ${T}_{{\rm{c}}}$, the correlation decays exponentially in time and the motion of the test star is diffusive (e.g., (Δ$\widehat{{\boldsymbol{L}}}$)² ∝ t). On these long timescales, the random walk of the test particle is analogous to that induced by fluctuations ${\delta }_{{\rm{D}}}$-correlated in time (as in the classical Brownian motion), leading to a diffusive random walk on the sphere.

Because it involves an integral over ${\boldsymbol{K}}$ that has to be computed for every value of t, Equation (33) remains difficult to implement. Let us now present a simpler toy model to generate a stochastic motion on the sphere that would share correlation properties similar to those of Equation (33). As such, we will assume that the stochastic motion of the test particle is generated by an effective dipole Gaussian noise, and follows the Langevin equation

$$\begin{eqnarray}&&\displaystyle \frac{d\widehat{{\boldsymbol{L}}}}{{dt}}={{\rm{\Gamma }}}_{{\rm{t}}}\,{\boldsymbol{\eta }}(t)\times \widehat{{\boldsymbol{L}}},\end{eqnarray} \tag{ 35 }$$

where the Gaussian noise ${\boldsymbol{\eta }}(t)$ is a 3D vector of zero mean, $\langle {\eta }_{{\text{}}i}(t)\rangle =0$, and follows $\langle {\eta }_{i}(t)\,{\eta }_{j}(t^{\prime} )\rangle ={\delta }_{{ij}}\,{e}^{-{[(t-t^{\prime} )/{T}_{{\rm{c}}}^{{\rm{t}}}]}^{2}}$. We will then choose the amplitude ${{\rm{\Gamma }}}_{{\rm{t}}}$ and coherence time ${T}_{{\rm{c}}}^{{\rm{t}}}$ by matching the short- and long-timescale behavior of the test particle's correlation function with the ballistic and diffusive regimes of the generic result from Equation (33), an approach already used in Hamers et al. (2018).

Following the same steps as in Equation (29), we may compute the correlation function of a test particle whose dynamics is imposed by Equation (35). It reads

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}(t) & = & \displaystyle \frac{{\delta }_{\alpha }^{\beta }}{4\pi }\,{e}^{-\tfrac{{A}_{{{\ell }}_{\alpha }}}{2}{{\rm{\Gamma }}}_{{\rm{t}}}^{2}{\left({T}_{{\rm{c}}}^{{\rm{t}}}\right)}^{2}\chi (t/{T}_{{\rm{c}}}^{{\rm{t}}})}\\ & \simeq & \ \displaystyle \frac{{\delta }_{\alpha }^{\beta }}{4\pi }\times \left\{\begin{array}{cc}\,{e}^{-\tfrac{{A}_{{{\ell }}_{\alpha }}}{2}{{\rm{\Gamma }}}_{{\rm{t}}}^{2}{t}^{2}} & t\ll {T}_{{\rm{c}}}^{{\rm{t}}},\\ {e}^{-\tfrac{\sqrt{\pi }{A}_{{{\ell }}_{\alpha }}}{2}{{\rm{\Gamma }}}_{{\rm{t}}}^{2}{T}_{{\rm{c}}}^{{\rm{t}}}t} & t\gg {T}_{{\rm{c}}}^{{\rm{t}}},\end{array}\right.\end{array}\end{eqnarray} \tag{ 36 }$$

with A_ℓ = ℓ(ℓ + 1). By matching the ballistic and diffusive regimes of Equation (36) with those of Equation (33), we may then constrain the amplitude, ${{\rm{\Gamma }}}_{{\rm{t}}}$, and coherence time, ${T}_{{\rm{c}}}^{{\rm{t}}}$, of the toy model of Equation (35).

The amplitude, ${{\rm{\Gamma }}}_{{\rm{t}}}$(${{\boldsymbol{K}}}_{{\rm{t}}}$), varies from realization to realization and is given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{t}}}^{2}({{\boldsymbol{K}}}_{{\rm{t}}})=\displaystyle \sum _{{\ell }}{B}_{{\ell }}\int d{\boldsymbol{K}}\,{n}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}})\,{{ \mathcal J }}_{{\ell }}^{2}\left[{{\boldsymbol{K}}}_{{\rm{t}}},{\boldsymbol{K}}\right].\end{eqnarray} \tag{ 37 }$$

When ensemble-averaged over realizations, this amplitude becomes

$$\begin{eqnarray}&&\langle {{\rm{\Gamma }}}_{{\rm{t}}}^{2}({{\boldsymbol{K}}}_{{\rm{t}}})\rangle ={{\rm{\Gamma }}}^{2}({{\boldsymbol{K}}}_{{\rm{t}}}),\end{eqnarray} \tag{ 38 }$$

as already defined in Equation (15). Finally, the coherence time is given by

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl} \langle {{\rm{\Gamma }}}_{{\rm{t}}}^{2}({{\boldsymbol{K}}}_{{\rm{t}}}) \rangle \,{T}_{{\rm{c}}}^{{\rm{t}}}({{\boldsymbol{K}}}_{{\rm{t}}}) & = & \displaystyle \sum _{\ell }\displaystyle \frac{{B}_{\ell }}{\sqrt{{A}_{\ell }/2}}\,\,\int d{\boldsymbol{K}}\,n({\boldsymbol{K}}){{J}}_{\ell }^{2}\left[{{\boldsymbol{K}}}_{{\rm{t}}},{\boldsymbol{K}}\right]\,{T}_{{\rm{c}}}({\boldsymbol{K}}),\end{array}\end{array}\,\end{eqnarray} \tag{ 39 }$$

where for simplicity we assume, similarly to Equation (23), that the coherence time, ${T}_{{\rm{c}}}^{{\rm{t}}}$(${{\boldsymbol{K}}}_{{\rm{t}}}$), is independent of the considered realization.

Equation (36) is the key result of this section. Indeed, it provides us with an analytical description of the statistical properties of the random walk of a test particle's orientation, as jointly induced by all the background particles. The test particle's random walk is characterized by the two quantities $({{\rm{\Gamma }}}_{{\rm{t}}}^{2},{T}_{{\rm{c}}}^{{\rm{t}}})$, which both depend on the test particle's parameters ${{\boldsymbol{K}}}_{{\rm{t}}}$. On the one hand, the torque amplitude, ${{\rm{\Gamma }}}_{{\rm{t}}}$, controls the amplitude of the test particle's initial ballistic motion. As such, the torque time, 1/${{\rm{\Gamma }}}_{{\rm{t}}}$, represents the typical time it would take the test particle to explore the unit sphere, for a given and frozen value of the background noise. On the other hand, the coherence time, ${T}_{{\rm{c}}}^{{\rm{t}}}$, describes the typical time one has to wait for the torque amplitude to reach a statistically different value, i.e., to decorrelate itself. It corresponds therefore to the timescale after which the test particle leaves the ballistic regime to enter the diffusive regime.

One strength of the present formalism is that, following Equations (37) and (39), one now has at one's disposal explicit expressions for these two parameters. These coefficients are time-independent, and can then easily be computed for various cluster models (by varying the DF n(${\boldsymbol{K}}$)) and various test particles (by varying ${{\boldsymbol{K}}}_{{\rm{t}}}$). Following Equation (36), the present approach also offers an explicit expression for the time-dependence of the test particle's correlation function, ${C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}(t)$. Considering the same test particle as in Figure 4, we illustrate in Figure 5 one random walk generated using the Langevin Equation (35).

Zoom In Zoom Out Reset image size

However, as we had already emphasized in Equation (22), it is important to note that the present system suffers from being non-ergodic, i.e., ensemble averages and time averages cannot be interchanged. This is highlighted by the fact that even for the exact same test particle (i.e., same ${{\boldsymbol{K}}}_{{\rm{t}}}$), the amplitude ${{\rm{\Gamma }}}_{{\rm{t}}}$ varies from realization to realization. In Appendix F, we compute the associated variance, as given by Equation (100), and show that this effect originates from the constraint of total energy conservation. Moreover, we show that this variance remains non-zero even in the limit of a Gaussian noise, and as such does not vanish in the limit of an infinite number of background particles generating the fluctuations.

Let us finally use our fiducial numerical simulations to highlight the result from Equation (36). This is illustrated in Figure 6 and, to shorten the main text, we detail in Appendix H.2 the procedure followed to obtain that figure. In that figure we note that the numerical measurements and the analytical prediction from Equation (36) agree both on short timescales but also on timescales longer than the coherence time ${T}_{{W}_{{\rm{t}}}}^{{\rm{t}}}$ (defined in Equation (118)). As already stressed in Equation (33), Figure 6 clearly exhibits the two successive regimes of evolution, namely ballistic for t ≪ ${T}_{{\rm{c}}}^{{\rm{t}}}$ and diffusive for t ≫ ${T}_{{\rm{c}}}^{{\rm{t}}}$. This same figure also emphasizes the importance of accounting for the variance in ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$, to correctly capture the late-time behavior of the test particles' stochastic motions. We recall that this effect that does not vanish in the limit of an infinite number of background particles. Since ${\varphi }_{(1,m)}^{{\rm{t}}}$(t) ∼ ${\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}$, Figure 6 also offers an illustration of the behavior of $t\mapsto \langle {\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}(t)\,{\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}(0)\rangle$. It is also straightforward to adapt that prediction to different test stars (by changing ${{\boldsymbol{K}}}_{{\rm{t}}}$) or different galactic nuclei (by changing the DF n(${\boldsymbol{K}}$)).

Zoom In Zoom Out Reset image size

In the previous derivations, we proceeded in two successive steps. First, in Section 3, we used estimates of the derivatives of the correlation function of the noise at the initial time to obtain an ansatz in Equation (18) for the time-dependence of the correlation function of the noise generated by the N background particles. Then, in Section 4, we used this noise as a source term to study the stochastic dynamics of a test particle. Yet, if the considered test particle is taken to be one particular background particle, its random walk in orientation and the background fluctuations sourcing it have to satisfy some self-consistency relation. This is what we explore in this section.

We start from Equation (20), and replace ${\varphi }_{\alpha }({\boldsymbol{K}},t)$ by its definition in terms of a discrete sum over particles, as in Equation (8), so that

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\alpha \beta }^{{\rm{r}}}({\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} },t-{t}^{{\prime} }) & = & \displaystyle \sum _{i,j}{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}_{i})\,{\delta }_{{\rm{D}}}({{\boldsymbol{K}}}^{{\prime} }-{{\boldsymbol{K}}}_{j})\\ & & \times \ {\left\langle {Y}_{\alpha }({\widehat{{\boldsymbol{L}}}}_{i}(t)){Y}_{\beta }({\widehat{{\boldsymbol{L}}}}_{j}({t}^{{\prime} }))\right\rangle }_{T}.\end{array}\end{eqnarray} \tag{ 40 }$$

We now assume that each background particle can be treated as a test particle, and that their long-term motions are decorrelated one from another. Only contributions from i = j remain, and Equation (40) becomes

$$\begin{eqnarray}&&{C}_{\alpha \beta }^{{\rm{r}}}({\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} },t)={\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}^{{\prime} })\displaystyle \sum _{i}{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}_{i})\,{C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}({{\boldsymbol{K}}}_{i},t),\end{eqnarray} \tag{ 41 }$$

where ${C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}({{\boldsymbol{K}}}_{i},t)$ is the test particle's correlation function of the particle i, as defined in Equation (31). To proceed further, let us now take the ensemble average of both sides of Equation (41), to get

$$\begin{eqnarray}&&{C}_{\alpha \beta }({\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} },t)={\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}^{{\prime} })\,4\pi \,n({\boldsymbol{K}})\,\left\langle {C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}({\boldsymbol{K}},t)\right\rangle .\end{eqnarray} \tag{ 42 }$$

Luckily, in Equation (32), we have already solved for the correlation function of the random walks of the test particle through the Magnus series. Using Equation (74), we generically obtain

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}({\boldsymbol{K}},t) & = & \displaystyle \frac{{\delta }_{\alpha }^{\beta }}{4\pi }\,\exp \left\{-\displaystyle \frac{{A}_{{{\ell }}_{\alpha }}}{2}\,\displaystyle \sum _{{\ell }}{B}_{{\ell }}\displaystyle \int d{{\boldsymbol{K}}}^{{\prime} }\,{{ \mathcal J }}_{{\ell }}^{2}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} }\right]\right.\\ & & \times \ \left.{\displaystyle \int }_{0}^{t}{{dt}}_{1}{\displaystyle \int }_{0}^{t}{{dt}}_{2}\,{C}_{{\ell }}^{{\rm{r}}}({{\boldsymbol{K}}}^{{\prime} },{t}_{1}-{t}_{2})\right\}.\end{array}\end{eqnarray} \tag{ 43 }$$

We may then take the ensemble-average of this relation and, for simplicity, keep only the first cumulant in the cumulant theorem. Reinjected into Equation (42), this leads to

$$\begin{eqnarray}\begin{array}{rcl}{C}_{{\ell }}({\boldsymbol{K}},t) & = & n({\boldsymbol{K}})\exp \left\{-\displaystyle \frac{{A}_{{\ell }}}{2}\displaystyle \sum _{{{\ell }}^{{\prime} }}{B}_{{{\ell }}^{{\prime} }}\displaystyle \int d{{\boldsymbol{K}}}^{{\prime} }\,{{ \mathcal J }}_{{{\ell }}^{{\prime} }}^{2}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} }\right]\right.\\ & & \times \ \left.{\displaystyle \int }_{0}^{t}{{dt}}_{1}{\displaystyle \int }_{0}^{t}{{dt}}_{2}\,{C}_{{{\ell }}^{{\prime} }}({{\boldsymbol{K}}}^{{\prime} },{t}_{1}-{t}_{2})\right\}.\end{array}\end{eqnarray} \tag{ 44 }$$

Equation (44) takes the form a self-consistent integral equation satisfied by the correlation of the noise fluctuations in the system. This relation can be further clarified by defining

$$\begin{eqnarray}&&{R}_{{\ell }}({\boldsymbol{K}},t)=\displaystyle \frac{1}{2}{\int }_{0}^{t}{{dt}}_{1}{\int }_{0}^{t}{{dt}}_{2}\,{C}_{{\ell }}({\boldsymbol{K}},{t}_{1}-{t}_{2}),\end{eqnarray} \tag{ 45 }$$

and one finally gets the self-consistent differential equation⁴

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{d}^{2}{R}_{{\ell }}({\boldsymbol{K}},t)}{{{dt}}^{2}} & = & n({\boldsymbol{K}})\,\exp \left\{-{A}_{{\ell }}\displaystyle \sum _{{{\ell }}^{{\prime} }}{B}_{{{\ell }}^{{\prime} }}\displaystyle \int d{{\boldsymbol{K}}}^{{\prime} }\,{{ \mathcal J }}_{{{\ell }}^{{\prime} }}^{2}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} }\right]\right.\\ & & \times \ \left.{R}_{{{\ell }}^{{\prime} }}({{\boldsymbol{K}}}^{{\prime} },t)\Space{0ex}{3.5ex}{0ex}\Space{0ex}{1.0ex}{0ex}\right\}.\end{array}\end{eqnarray} \tag{ 46 }$$

Equation (46) is the important result of this section, as it highlights the self-consistent relation satisfied by the correlation of the noise fluctuations. Yet, as it couples both different harmonics (via ${\sum }_{{{\ell }}^{{\prime} }}$) and different parameters (via $\int d{{\boldsymbol{K}}}^{{\prime} }$), such a differential equation appears too intricate to easily be solved explicitly. This is not pursued further here.

One may still proceed iteratively to obtain improved approximations of the noise correlation function. To do so, one starts from the Gaussian dependence obtained in Equation (18). This (motivated) ansatz may then be reinjected in the rhs of the self-consistency relation from Equation (44), leading to a new expression of the noise correlation function that would have both a ballistic and a diffusive part. Such a procedure is illustrated in Figure 3, where we show how one can better match the late-time properties of the system's noise through this iterative process.

As an illustration of the present formalism, let us now consider the case of a stellar cusp distribution similar to that of SgrA*. The mass of the MBH is taken to be ${M}_{\bullet }$ = 4 × 10⁶ M_⊙, and for simplicity we consider a single-mass stellar population of individual mass m_⋆ = 1 M_⊙. We assume that the stars' eccentricities follow a thermal distribution, f_e(e) = 2e (Merritt 2013), and that the number of stars per unit a follows a power-law distribution of the form ${n}_{a}(a)=({N}_{0}/{a}_{0}){(a/{a}_{0})}^{2-\gamma }$, where N₀ = g(γ)N(<a₀), with

$$\begin{eqnarray}&&g(\gamma )={2}^{-\gamma }\,(3-\gamma )\,\sqrt{\pi }\,\displaystyle \frac{{\rm{\Gamma }}(1+\gamma )}{{\rm{\Gamma }}\left(\gamma -{\textstyle \tfrac{1}{2}}\right)},\end{eqnarray} \tag{ 47 }$$

and N(<a₀) the physical number of stars within a sphere of radius a₀ from the center. For the numerical application, we assume that a₀ = ${r}_{{\rm{h}}}$ = 2 $\mathrm{pc}$ and N(<a₀) = 4 × 10⁶. The system being of infinite extent, we write the system's DF as n(m, a, e) = f_m(m)f_e(e) n_a(a)/(4π).

In Appendix I, we show that the amplitude Γ², defined in Equation (15) and characterizing the ballistic regime of the orientation's random walk, follows the power-law distribution

$$\begin{eqnarray}&&{{\rm{\Gamma }}}^{2}(a,e)=\displaystyle \frac{N(\lt a)}{{P}^{2}(a)}\displaystyle \frac{\langle {m}^{2}\rangle }{{M}_{\bullet }^{2}}\displaystyle \frac{\pi g(\gamma ){f}_{{{\rm{\Gamma }}}^{2}}(e)}{1-{e}^{2}},\end{eqnarray} \tag{ 48 }$$

where P(a) = 2π(a³/${({{GM}}_{\bullet }))}^{1/2}$ is the orbital period, and ${f}_{{{\rm{\Gamma }}}^{2}}(e)\simeq 0.15$ is a dimensionless eccentricity function defined in Equation (127) and illustrated in Figure 10. In Figure 7, we illustrate the dependence of the torque time, 1/Γ, for circular orbits of different semimajor axes and for different cusp profiles, and interestingly note that this VRR timescale is similar to the age of some of the young stars observed in our Galactic center (Habibi et al. 2017).

Zoom In Zoom Out Reset image size

As shown in Figure 7, should the S-stars be born in a disk, the VRR process is sufficiently fast to isotropize their orbital orientations (Hopman & Alexander 2006), but SRR may not be efficient enough to thermalize their eccentricities (Bar-Or & Fouvry 2018).

One can follow a similar calculation to obtain the expression of the coherence time, ${T}_{{\rm{c}}}^{{\rm{t}}}$, defined in Equation (39) and characterizing the diffusive regime of the orientation's random walk. It follows the power-law distribution

$$\begin{eqnarray}&&{T}_{{\rm{c}}}^{{\rm{t}}}(a,e)=\displaystyle \frac{P(a)}{\sqrt{N(\lt a)}}\displaystyle \frac{{M}_{\bullet }}{\sqrt{\langle {m}^{2}\rangle }}\,\displaystyle \frac{{f}_{T}(e)}{\sqrt{\pi g(\gamma )}},\end{eqnarray} \tag{ 49 }$$

where the dimensionless eccentricity function f_T (e) is given in Equation (132) and illustrated in Figure 10. In that figure, we note that for a thermal eccentricity distribution and a cusp's power index 1 ≤ γ ≤ 2, one can assume that ${f}_{{{\rm{\Gamma }}}^{2}}^{1/2}(e)\,{f}_{T}(e)\simeq 0.4$, which leads to the torque time and the coherence time following the approximate relation

$$\begin{eqnarray}&&{T}_{{\rm{c}}}^{{\rm{t}}}(a,e)\simeq \displaystyle \frac{0.4}{\sqrt{1-{e}^{2}}}\,\displaystyle \frac{1}{{\rm{\Gamma }}(a,e)}.\end{eqnarray} \tag{ 50 }$$

We note that this simple relation allows for an even simpler generation of samples of random walks in orientations as given by the toy model from Equation (35), as one only has to estimate the test particle's torque time 1/Γ (a, e), as the associated coherence time, ${T}_{{\rm{c}}}^{{\rm{t}}}$(a, e), follows immediately.

In the present work, we illustrated how one can describe quantitatively the statistical properties of the stochastic evolutions of star's orientations in galactic nuclei during the process of VRR. The main difficulty in the present derivation lies in the system being fundamentally degenerate, i.e., having a vanishing mean field Hamiltonian, H = 0. This system is also non-Markovian, i.e., correlated in time, as well as non-ergodic, i.e., time- and ensemble-averages cannot be interchanged.

Placing ourselves in the limit of an isotropic distribution of stars, we circumvented some of these difficulties in Section 3 by assuming that the statistical properties of the noise fluctuations can be derived from estimates of the derivatives of their correlation function at the initial time. The main result was obtained in Equation (22), which provided us with a self-consistent ansatz for the statistical properties of the time-dependence of the correlation of the fluctuations generated jointly by the system's N particles. In Section 4, we used this result to describe the random walk of a test particle's orientation embedded in this stochastic system, recovering both the ballistic and diffusive regimes. The main result was obtained in Equation (36), which yields quantitative predictions for the statistical properties of that random walk. The key tools used at that stage were the Magnus series to solve the linear matrix evolution equation for the test particle, the independence hypothesis to separate the statistics of the background noise from that of the test particle's random walk, and the cumulant theorem to estimate ensemble averages. We also emphasized how non-ergodic effects (associated with the constraint of total energy conservation) should be accounted for to allow for reliable long-timescale predictions. Throughout the text, all the predictions were compared with detailed numerical simulations offering a quantitative agreement. In Section 5, we highlighted the self-consistency existing between the spontaneous fluctuations in the system and the associated random walks in orientations. Finally, in Section 6, we presented a first application of this framework to estimate the timescales of VRR in a stellar cusp similar to that of SgrA*.

This paper is only a first step toward a complete theory of VRR, and we list below some possible avenues for future development. In the current derivation, we relied extensively on the isotropic assumption, and as such neglected any effects associated with anisotropic clustering in orientation (Szölgyén & Kocsis 2018). For binaries, the exact statistical properties of the VRR random walk in orientation can lead to enhanced rates of mergers (Hamers et al. 2018), hence the importance of quantitative predictions for the properties of these random walks, as obtained in Equation (36). Building upon Section 4, one could also investigate how a substructure like a disk stochastically dissolves (Kocsis & Tremaine 2011). This asks for a detailed accounting of the correlations in the potential fluctuations of a given realization, to characterize how stars with similar initial orientations or similar parameters get slowly separated. Finally, here we focused our interest on systems dominated by a central mass. Provided one updates accordingly the interaction coupling coefficients, ${{ \mathcal J }}_{{\ell }}$[${\boldsymbol{K}}$, ${{\boldsymbol{K}}}^{{\prime} }$], similar investigations could be pursued in the context of spherical globular clusters (Meiron & Kocsis 2018).

We thank Christophe Pichon and Scott Tremaine for remarks on an earlier version of this manuscript. J.B.F. acknowledges support from Program number HST-HF2-51374 which was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5–26555. B.B. is supported by membership from Martin A. and Helen Chooljian at the Institute for Advanced Study.

Similarly to Equation (10) in Kocsis & Tremaine (2015), we define the coupling coefficients⁵

$$\begin{eqnarray}&&{{ \mathcal J }}_{{\ell }}\left[{{\boldsymbol{K}}}_{i},{{\boldsymbol{K}}}_{j}\right]=\displaystyle \frac{{{Gm}}_{i}{m}_{j}}{{a}_{\mathrm{out}}}\displaystyle \frac{1}{L[{{\boldsymbol{K}}}_{i}]}\,{s}_{{\ell }}\left[\alpha ,{e}_{\mathrm{in}},{e}_{\mathrm{out}}\right],\end{eqnarray} \tag{ 51 }$$

where $L[{\boldsymbol{K}}]=m\sqrt{{{GM}}_{\bullet }a(1-{e}^{2})}$ is the magnitude of the angular momentum. We also introduce "out" (resp. "in") as the index i or j with the larger (resp. smaller) semimajor axis, and define accordingly the ratio α = ${a}_{\mathrm{in}}$/${a}_{\mathrm{out}}$ ≤ 1. With this notation, the dimensionless coefficients s_ℓ[α, ${e}_{\mathrm{in}}$, ${e}_{\mathrm{out}}$] are given by

$$\begin{eqnarray}&&\begin{array}{l}{s}_{{\ell }}\left[\alpha ,{e}_{\mathrm{in}},{e}_{\mathrm{out}}\right]=\displaystyle \frac{4\pi \,{P}_{{\ell }}^{2}(0)}{2{\ell }+1}\,\displaystyle \frac{1}{\alpha }\,\displaystyle \frac{1}{{\pi }^{2}}{\displaystyle \int }_{0}^{\pi }d\phi {\displaystyle \int }_{0}^{\pi }d{\phi }^{{\prime} }\\ \times \ \displaystyle \frac{{\left(\min \left[\alpha (1-{e}_{\mathrm{in}}\cos (\phi )),(1-{e}_{\mathrm{out}}\cos ({\phi }^{{\prime} }))\right]\right)}^{{\ell }+1}}{{\left(\max \left[\alpha (1-{e}_{\mathrm{in}}\cos (\phi )),(1-{e}_{\mathrm{out}}\cos ({\phi }^{{\prime} }))\right]\right)}^{{\ell }}},\end{array}\end{eqnarray} \tag{ 52 }$$

with P_ℓ(u) the usual Legendre polynomials. Because they are independent of the details of the considered system, the coefficients s_ℓ[α, ${e}_{\mathrm{in}}$, ${e}_{\mathrm{out}}$] can be precomputed on a grid to hasten the numerical evaluation of ${{ \mathcal J }}_{{\ell }}$[K_i, K_j]. For our fiducial simulations, these coefficients were pre-computed on a linear 3D grid in (α, ${e}_{\mathrm{in}}$, ${e}_{\mathrm{out}}$) consisting of 200³ elements, with 10⁻² ≤ α ≤ 1 and 0 ≤ ${e}_{\mathrm{in}}$, ${e}_{\mathrm{out}}$ ≤ 0.99. We refer to Figure 1 in Kocsis & Tremaine (2015) for an illustration of the behavior of these coefficients.

In this appendix, we follow James (1973) and Ivers & Phillips (2008) and detail some of the properties of the Elsässer coefficients. We emphasize that we work with real spherical harmonics, hence the need for some identities to be adapted.

The real spherical harmonics are defined with the convention

$$\begin{eqnarray}&&\begin{array}{l}{Y}_{{\ell }m}(u,\phi )\\ \ =\ \left\{\begin{array}{cc}\sqrt{2}\,{K}_{{\ell }}^{| m| }\,{P}_{{\ell }}^{| m| }(u)\,\sin (| m| \phi ) & \mathrm{if}\,m\lt 0,\\ {K}_{{\ell }}^{0}\,{P}_{{\ell }}^{0}(u) & \mathrm{if}\,m=0,\\ \sqrt{2}\,{K}_{{\ell }}^{m}\,{P}_{{\ell }}^{m}(u)\,\cos (m\phi ) & \mathrm{if}\,m\gt 0,\end{array}\right.\end{array}\end{eqnarray} \tag{ 53 }$$

with ${P}_{{\ell }}^{m}$(u) the usual associated Legendre polynomials, and the coefficients ${K}_{{\ell }}^{m}={\left[\tfrac{2{\ell }+1}{4\pi }\tfrac{({\ell }-m)!}{({\ell }+m)!}\right]}^{1/2}$. With this convention, the spherical harmonics follow the normalization $\int d\widehat{{\boldsymbol{L}}}\,{Y}_{{\ell }m}{Y}_{{{\ell }}^{{\prime} }{m}^{{\prime} }}$ = ${\delta }_{{\ell }}^{{{\ell }}^{{\prime} }}{\delta }_{m}^{{m}^{{\prime} }}$.

Following Equation (5.5) of Ivers & Phillips (2008), the real Elsässer coefficients, as defined in Equation (10), can be decomposed as

$$\begin{eqnarray}&&{E}_{\alpha \gamma \delta }={E}_{{{\ell }}_{\alpha }{{\ell }}_{\gamma }{{\ell }}_{\delta }}^{L}\,{E}_{\alpha \gamma \delta }^{M},\end{eqnarray} \tag{ 54 }$$

where ${E}_{{{\ell }}_{\alpha }{{\ell }}_{\gamma }{{\ell }}_{\delta }}^{L}$ only depends on (ℓ_α, ℓ_γ, ℓ_δ), while ${E}_{\alpha \gamma \delta }^{M}$ also depends on (m_α, m_γ, m_δ). The m-independent coefficients read

$$\begin{eqnarray}{E}_{{{\ell }}_{\alpha }{{\ell }}_{\gamma }{{\ell }}_{\delta }}^{L}={{\rm{\Lambda }}}_{{{\ell }}_{\alpha }{{\ell }}_{\gamma }{{\ell }}_{\delta }}\,{{\rm{\Delta }}}_{{{\ell }}_{\alpha }{{\ell }}_{\gamma }{{\ell }}_{\delta }}\left(\begin{array}{ccc}{{\ell }}_{\alpha }+1 & {{\ell }}_{\gamma }+1 & {{\ell }}_{\delta }+1\\ 0 & 0 & 0\end{array}\right),\end{eqnarray} \tag{ 55 }$$

where we introduce the Wigner 3j-symbols (Arfken et al. 2005), and define

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Lambda }}}_{{{\ell }}_{\alpha }{{\ell }}_{\gamma }{{\ell }}_{\delta }}=\sqrt{\displaystyle \frac{(2{{\ell }}_{\alpha }+1)(2{{\ell }}_{\gamma }+1)(2{{\ell }}_{\delta }+1)}{4\pi }},\\ {{\rm{\Delta }}}_{{{\ell }}_{\alpha }{{\ell }}_{\gamma }{{\ell }}_{\delta }}=\sqrt{\displaystyle \frac{({{\ell }}_{\alpha }+{{\ell }}_{\gamma }+{{\ell }}_{\delta }+2)({{\ell }}_{\alpha }+{{\ell }}_{\gamma }+{{\ell }}_{\delta }+4)}{4({{\ell }}_{\alpha }+{{\ell }}_{\gamma }+{{\ell }}_{\delta }+3)}}\\ \times \ \sqrt{({{\ell }}_{\alpha }+{{\ell }}_{\gamma }-{{\ell }}_{\delta }+1)({{\ell }}_{\delta }+{{\ell }}_{\alpha }-{{\ell }}_{\gamma }+1)({{\ell }}_{\gamma }+{{\ell }}_{\delta }-{{\ell }}_{\alpha }+1)}.\end{array}\end{eqnarray} \tag{ 56 }$$

The coefficients ${E}_{\alpha \gamma \delta }^{M}$ are given by

$$\begin{eqnarray}{E}_{\alpha \gamma \delta }^{M}=\displaystyle \sum _{{\varepsilon }_{\alpha },{\varepsilon }_{\gamma },{\varepsilon }_{\delta }=\pm 1}\,\mathrm{Im}\left[{K}_{\alpha \gamma \delta }^{{\varepsilon }_{\alpha }{\varepsilon }_{\gamma }{\varepsilon }_{\delta }}\right]\,\left(\begin{array}{ccc}{{\ell }}_{\alpha } & {{\ell }}_{\gamma } & {{\ell }}_{\delta }\\ {\varepsilon }_{\alpha }{m}_{\alpha } & {\varepsilon }_{\gamma }{m}_{\gamma } & {\varepsilon }_{\delta }{m}_{\delta }\end{array}\right),\end{eqnarray} \tag{ 57 }$$

where the tensor ${K}_{\alpha \gamma \delta }^{{\varepsilon }_{\alpha }{\varepsilon }_{\gamma }{\varepsilon }_{\delta }}$ comes from the fact that we are considering real spherical harmonics, and is given by

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{K}_{\alpha \gamma \delta }^{{\varepsilon }_{\alpha }{\varepsilon }_{\gamma }{\varepsilon }_{\delta }} & = & {\kappa }_{{m}_{\alpha }}^{{\varepsilon }_{\alpha }}\,{\kappa }_{{m}_{\gamma }}^{{\varepsilon }_{\gamma }}\,{\kappa }_{{m}_{\delta }}^{{\varepsilon }_{\delta }}\,\,{\rm{w}}{\rm{i}}{\rm{t}}{\rm{h}}\\ {\kappa }_{m}^{+1} & = & \left\{\begin{array}{cc}i\,{\left(-1\right)}^{m}/\sqrt{2} & {\rm{i}}{\rm{f}}\,m\lt 0\\ 1/2 & {\rm{i}}{\rm{f}}\,m=0\\ 1/\sqrt{2} & {\rm{i}}{\rm{f}}\,m\gt 0\end{array}\right.\,{\rm{a}}{\rm{n}}{\rm{d}}\\ {\kappa }_{m}^{-1} & = & \left\{\begin{array}{cc}-i/\sqrt{2} & {\rm{i}}{\rm{f}}\,m\lt 0\\ 1/2 & {\rm{i}}{\rm{f}}\,m=0\\ {\left(-1\right)}^{m}/\sqrt{2} & {\rm{i}}{\rm{f}}\,m\gt 0.\end{array}\right.\end{array}\end{array}\end{eqnarray} \tag{ 58 }$$

The Elsässer coefficients satisfy various exclusion rules (James 1973). In particular, for E_αγδ to be non-zero, one has to satisfy

$$\begin{eqnarray}&&\begin{array}{l}\{{\rm{C}}1\}:\,| {m}_{\alpha }| \leqslant {{\ell }}_{\alpha };\,| {m}_{\gamma }| \leqslant {{\ell }}_{\gamma };\,| {m}_{\delta }| \leqslant {{\ell }}_{\delta },\\ \{{\rm{C}}2\}:\,{{\ell }}_{\alpha }+{{\ell }}_{\gamma }+{{\ell }}_{\delta }\ \mathrm{is}\ \mathrm{odd},\\ \{{\rm{C}}3\}:\,| {{\ell }}_{\alpha }-{{\ell }}_{\gamma }| \lt {{\ell }}_{\delta }\lt {{\ell }}_{\alpha }+{{\ell }}_{\gamma }\ (\mathrm{strict}\ \mathrm{triangular}\ \mathrm{inequality}),\\ \{{\rm{C}}4\}:\,\mathrm{all}\ \mathrm{pairs}\ ({{\ell }}_{\alpha },{m}_{\alpha }),({{\ell }}_{\gamma },{m}_{\gamma }),({{\ell }}_{\delta },{m}_{\delta })\ \mathrm{are}\ \mathrm{different}.\end{array}\end{eqnarray} \tag{ 59 }$$

These coefficients also follow the symmetry relations E_αδγ = E_γαδ = −E_αγδ.

Finally, following Varshalovich et al. (1988), the Elsässer coefficients satisfy various contraction identities. In particular, throughout the derivations, we will rely on

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{{m}_{\gamma },{m}_{\delta }}\,{E}_{\alpha \gamma \delta }^{M}\,{E}_{\beta \gamma \delta }^{M}={\delta }_{\alpha }^{\beta }\,\displaystyle \frac{1}{2{{\ell }}_{\alpha }+1}\ \mathrm{for}\ | {{\ell }}_{\alpha }-{{\ell }}_{\beta }| \lt {{\ell }}_{\gamma /\delta }\lt | {{\ell }}_{\alpha }+{{\ell }}_{\beta }| \ \mathrm{and}\ {{\ell }}_{\alpha }+{{\ell }}_{\beta }+{{\ell }}_{\gamma /\delta }\,\mathrm{odd},\\ \displaystyle \sum _{{{\ell }}_{\delta }}\displaystyle \frac{1}{2{{\ell }}_{\alpha }+1}{\left({E}_{{{\ell }}_{\alpha }{{\ell }}_{\gamma }{{\ell }}_{\delta }}^{L}\right)}^{2}={A}_{{{\ell }}_{\alpha }}\,{B}_{{{\ell }}_{\gamma }}\,\ \mathrm{with}\,{A}_{{\ell }}\,=\,{\ell }({\ell }+1),\,{B}_{{\ell }}=\displaystyle \frac{{\ell }({\ell }+1)(2{\ell }+1)}{8\pi }.\end{array}\end{eqnarray} \tag{ 60 }$$

In this appendix, we briefly detail the numerical simulations to which our analytical results are compared. To simulate a system of N interacting particles, the starting point is the evolution Equation (4), that can be used for each of the N particles. In that form, we note that the velocity vector, $d{\widehat{{\boldsymbol{L}}}}_{i}$/$d$t, is expressed only as a function of the current location of the particle, $\widehat{{\boldsymbol{L}}}$_i(t), and the instantaneous particle's magnetizations, M_ℓm (K_i, t). There are N such evolution equations, but because the magnetizations vary from one particle to another, their computation has to be made once per timestep and particle. As a result, the overall complexity of advancing the particles for one timestep scales like ${ \mathcal O }({N}^{2}{{\ell }}_{\max }^{2})$, with ℓ_max the maximum harmonic number considered in the pairwise interaction.⁶ In our approach, the motion of the particles is integrated by computing the magnetizations, while in the implementation presented in Kocsis & Tremaine (2015), particles are moved forward by solving successively pairwise interactions, an approach symplectic by design. The method of Kocsis & Tremaine benefits from parallelization by computing the pairwise interactions in parallel, provided that they are carefully ordered; see their Section 3.2. The present approach benefits from parallelization when computing the ${ \mathcal O }(N{{\ell }}_{\max }^{2})$ magnetizations, seen as a contraction of large matrices and vectors.

Our numerical implementation proceeds then by (i) computing efficiently the spherical harmonics (and the vector ones) at the location of the particles, (ii) computing the magnetizations in Equation (2), (iii) computing the velocity fields in Equation (4), and (iv) advancing all the particles' orientation for one timestep. The real spherical harmonics are computed using a recurrence relation for the renormalized associated Legendre polynomials (see Equation (6.7.9) in Press et al. 2007), and using the second-order recurrence relation $\cos (m\phi )\,=2\cos (\phi )\cos ((m-1)\phi )-\cos ((m-2)\phi )$ (similarly for sin(m ϕ)) for the azimuthal component. To compute the real vector spherical harmonics, we follow the same recurrence as in Appendix B.2 of Mignard & Klioner (2012), adapted to the renormalized associated Legendre polynomials. Once all the velocity vectors $d{\widehat{{\boldsymbol{L}}}}_{i}$/$d$t are computed, particles are advanced for a timestep h, using a fourth-order Runge–Kutta integrator (see Equation (17.1.3) in Press et al. 2007).

All the derivations presented in the main text are illustrated with comparisons using this numerical approach. We consider a system composed of N = 10³ stars, and assume that the particles' conserved quantities K_i = (m_i, a_i, e_i) satisfy m = ${m}_{\min }$, ${a}_{\min }$ ≤ a ≤ ${a}_{\max }$, and ${e}_{\min }$ ≤ e ≤ ${e}_{\max }$. Our units are chosen so that ${m}_{\min }$ = ${a}_{\min }$ = G = 1, and we pick ${a}_{\max }$/${a}_{\min }$ = 100, ${e}_{\min }$ = 0, and ${e}_{\max }$ = 0.3. These parameters are drawn independently one from another, according to probability distribution functions (PDFs) proportional to (${\delta }_{{\rm{D}}}$(m − ${m}_{\min }$), a^1/2, e), which corresponds to a single-mass population in a harmonic profile with a thermal distribution of small eccentricities. The stars' initial orientations are drawn uniformly on the sphere, and the interactions are truncated at ℓ_max = 50. The timestep of the simulation is the same for all particles and is determined at the start of each realization. To do so, we compute the torque exerted on every particle at the initial time, ${\tau }_{i}=| d{\widehat{{\boldsymbol{L}}}}_{i}/{dt}| =1/{t}_{\tau }^{i}$, and define ${t}_{\tau }^{i}$ as the associated torque time. The integration timestep is then fixed initially to $h={10}^{-2}\times {\mathrm{Min}}_{i}[{t}_{\tau }^{i}]$. With such choices, integrating the system for one timestep takes approximately 1 s on a single core, and simulations are carried out for 2 × 10⁵ timesteps.

D.1. Computing Ensemble Averages

We are generically interested in computing ensemble averages at the initial time of the form $\langle {\varphi }_{\alpha }({\boldsymbol{K}},0)\,{\varphi }_{\beta }({{\boldsymbol{K}}}^{{\prime} },0)...\rangle$. Such averages can be carried out explicitly by noting that, at the initial time, the N particles are drawn independently one from another, both for their orientations and their parameters. Following our isotropic assumption, their orientation is drawn uniformly on the sphere, according to the PDF f ($\widehat{{\boldsymbol{L}}}$) = 1/(4π), while we assume that their parameter ${\boldsymbol{K}}$ is drawn according to a PDF g(${\boldsymbol{K}}$), normalized so that $\int d{\boldsymbol{K}}\,g({\boldsymbol{K}})=1$.

To illustrate the gist of these calculations, let us consider the case $\langle {\varphi }_{\alpha }({\boldsymbol{K}},0)\,{\varphi }_{\beta }({{\boldsymbol{K}}}^{{\prime} },0)\rangle$. As they do not contribute to the dynamics, we never need to consider the harmonics (ℓ, m) = (0, 0), so that $\int d\widehat{{\boldsymbol{L}}}f(\widehat{{\boldsymbol{L}}}){Y}_{\alpha /\beta }(\widehat{{\boldsymbol{L}}})=0$. Owing to the particle independence at the initial time and following the definition from Equation (8), we can write

$$\begin{eqnarray}&&\begin{array}{l}\left\langle {\varphi }_{\alpha }({\boldsymbol{K}},0)\,{\varphi }_{\beta }({{\boldsymbol{K}}}^{{\prime} },0)\right\rangle \\ \ =\ \displaystyle \int d{\widehat{{\boldsymbol{L}}}}_{1}d{{\boldsymbol{K}}}_{1}...d{\widehat{{\boldsymbol{L}}}}_{N}d{{\boldsymbol{K}}}_{N}\,f({\widehat{{\boldsymbol{L}}}}_{1})\,g({{\boldsymbol{K}}}_{1})...f({\widehat{{\boldsymbol{L}}}}_{N})\,g({{\boldsymbol{K}}}_{N})\\ \ \ \ \times \ \displaystyle \sum _{i,j}{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}_{i})\,{\delta }_{{\rm{D}}}({{\boldsymbol{K}}}^{{\prime} }-{{\boldsymbol{K}}}_{j})\,{Y}_{\alpha }({\widehat{{\boldsymbol{L}}}}_{i})\,{Y}_{\beta }({\widehat{{\boldsymbol{L}}}}_{j})\\ \ =\ N{\left\langle {\varphi }_{\alpha }({\boldsymbol{K}},0){\varphi }_{\beta }({{\boldsymbol{K}}}^{{\prime} },0)\right\rangle }_{{\rm{c}}}\end{array}\end{eqnarray} \tag{ 61 }$$

where non-zero terms only come from i = j, and we introduce the connected average as

$$\begin{eqnarray}\begin{array}{rcl}{\left\langle {\varphi }_{\alpha }({\boldsymbol{K}},0){\varphi }_{\beta }({{\boldsymbol{K}}}^{{\prime} },0)\right\rangle }_{{\rm{c}}} & = & \displaystyle \int d{\widehat{{\boldsymbol{L}}}}_{1}d{{\boldsymbol{K}}}_{1}\,f({\widehat{{\boldsymbol{L}}}}_{1})\,g({{\boldsymbol{K}}}_{1}){\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}_{1})\\ & & \times \ {\delta }_{{\rm{D}}}({{\boldsymbol{K}}}^{{\prime} }-{{\boldsymbol{K}}}_{1})\,{Y}_{\alpha }({\widehat{{\boldsymbol{L}}}}_{1})\,{Y}_{\beta }({\widehat{{\boldsymbol{L}}}}_{1})\\ & = & \ {\delta }_{\alpha }^{\beta }\,{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{{\boldsymbol{K}}}^{{\prime} })\,\displaystyle \frac{g({\boldsymbol{K}})}{4\pi }.\end{array}\end{eqnarray} \tag{ 62 }$$

When considering averages at the initial time involving more than two fields, we limit ourselves to the dominant contributions associated with pair couplings (i.e., the limit of Gaussian fields, Wick's theorem). We can then write

$$\begin{eqnarray}&&\begin{array}{l}\left\langle {\varphi }_{1}({{\boldsymbol{K}}}_{1},0)...{\varphi }_{n}({{\boldsymbol{K}}}_{n},0)\right\rangle \\ \ =\ \left\{\begin{array}{cc}{N}^{n/2}\left({\left\langle {\varphi }_{1}{\varphi }_{2}\right\rangle }_{{\rm{c}}}...{\left\langle {\varphi }_{n-1}{\varphi }_{n}\right\rangle }_{{\rm{c}}}+\mathrm{perm}.\right) & \ \mathrm{if}\ {\text{}}n\ \mathrm{even},\\ 0 & \ \mathrm{if}\ {\text{}}n\ \mathrm{odd},\end{array}\right.\end{array}\end{eqnarray} \tag{ 63 }$$

where "perm" browses all the possible pair decompositions without repetitions, and averages involving an odd number of fields are neglected.

D.2. Initial Values of the Correlation Function

Following the method just described, we may now estimate the value and the second derivative of the correlation function at the initial time, as introduced in Equation (11).

For the value at the initial time, we can write

$$\begin{eqnarray}\begin{array}{rcl}{C}_{\alpha \beta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,0) & = & \left\langle {\varphi }_{\alpha }({\boldsymbol{K}},0)\,{\varphi }_{\beta }({\boldsymbol{K}}^{\prime} ,0)\right\rangle \\ & = & \ {\delta }_{\alpha }^{\beta }\,{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{\boldsymbol{K}}^{\prime} )\,n({\boldsymbol{K}}),\end{array}\end{eqnarray} \tag{ 64 }$$

where we follow Equation (61) to compute the last average, and introduced n(${\boldsymbol{K}}$) = g(${\boldsymbol{K}}$) N/(4π) as the DF of the stars' parameters satisfying the normalization $\int d\widehat{{\boldsymbol{L}}}d{\boldsymbol{K}}\,n({\boldsymbol{K}})=N$. We emphasize that to compute Equation (64), we relied on the assumption of an isotropic distribution of particles on the sphere, which led to the Kronecker coefficients ${\delta }_{\alpha }^{\beta }$ w.r.t. the harmonic coefficients.

The ensemble average expectation for the first derivative at the initial time reads $\partial {C}_{\alpha \beta }/\partial t$ = $\langle {\dot{\varphi }}_{\alpha }\,{\varphi }_{\beta }\rangle$ ∼ $\langle {\varphi }_{\gamma }\,{\varphi }_{\delta }{\varphi }_{\beta }\rangle =0$, where we used the quadratic evolution Equation (9) once. As it involves an odd number of fields, this correlation is equal to zero, as imposed by Equation (63).

Let us now turn to the computation of the ensemble-average expectation for the second derivative of the correlation at the initial time. We write

$$\begin{eqnarray}&&\begin{array}{l}{\left.\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}{C}_{\alpha \beta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,t)\right|}_{t=0}\\ \ =\ \left\langle {\ddot{\varphi }}_{\alpha }({\boldsymbol{K}},0)\,{\varphi }_{\beta }({\boldsymbol{K}}^{\prime} ,0)\right\rangle \\ \ =\ -\left\langle {\dot{\varphi }}_{\alpha }({\boldsymbol{K}},0)\,{\dot{\varphi }}_{\beta }({\boldsymbol{K}}^{\prime} ,0)\right\rangle \\ \ =\ -\displaystyle \int d{{\boldsymbol{K}}}_{1}d{{\boldsymbol{K}}}_{2}\,{{ \mathcal J }}_{{{\ell }}_{\gamma }}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}_{1}\right]{{ \mathcal J }}_{{{\ell }}_{\delta }}\left[{\boldsymbol{K}}^{\prime} ,{{\boldsymbol{K}}}_{2}\right]{E}_{\alpha \gamma \lambda }{E}_{\beta \delta \rho }\\ \ \ \ \times \ \left\langle {\varphi }_{\gamma }({{\boldsymbol{K}}}_{1},0)\,{\varphi }_{\lambda }({\boldsymbol{K}},0)\,{\varphi }_{\delta }({{\boldsymbol{K}}}_{2},0)\,{\varphi }_{\rho }({\boldsymbol{K}}^{\prime} ,0)\right\rangle \end{array}\end{eqnarray} \tag{ 65 }$$

where we inject the evolution Equation (9) twice. As shown in Equation (63), in the limit of Gaussian fluctuations, the average term can be computed by keeping only averages of pairs. As E_αγγ = 0, only two of the possible couplings remain, namely $\left\langle \gamma \delta \right\rangle$ $\left\langle \lambda \rho \right\rangle$ and $\left\langle \gamma \rho \right\rangle$ $\left\langle \delta \lambda \right\rangle$, which leads to

$$\begin{eqnarray}&&\begin{array}{l}{\left.\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}{C}_{\alpha \beta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,t)\right|}_{t=0}\\ \quad \ =\ -{E}_{\alpha \gamma \delta }\,{E}_{\beta \gamma \delta }\left[{{\rm{\Gamma }}}_{\gamma }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} )-{{\rm{\Lambda }}}_{\gamma \delta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} )\right],\end{array}\end{eqnarray} \tag{ 66 }$$

where we use E_βδγ = −E_βγδ, and introduce

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}_{\gamma }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ) & = & {\delta }_{{\rm{D}}}({\boldsymbol{K}}-{\boldsymbol{K}}^{\prime} )\,n({\boldsymbol{K}})\displaystyle \int d{{\boldsymbol{K}}}^{{\prime\prime} }\,n({{\boldsymbol{K}}}^{{\prime\prime} })\,{{ \mathcal J }}_{{{\ell }}_{\gamma }}^{2}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime\prime} }\right],\\ {{\rm{\Lambda }}}_{\gamma \delta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ) & = & n({\boldsymbol{K}})\,n({\boldsymbol{K}}^{\prime} )\,{{ \mathcal J }}_{{{\ell }}_{\gamma }}\left[{\boldsymbol{K}},{\boldsymbol{K}}^{\prime} \right]\,{{ \mathcal J }}_{{{\ell }}_{\delta }}\left[{\boldsymbol{K}}^{\prime} ,{\boldsymbol{K}}\right].\end{array}\end{eqnarray} \tag{ 67 }$$

Following Equation (60), one can now perform the sums over m_γ and m_δ in Equation (66). The VRR interactions being limited to even harmonic numbers ℓ, we may impose at this stage that ℓ_α is even. Glancing back at the constraint {C2} from Equation (59), we note that ℓ_α + ℓ_γ + ℓ_δ has to be odd, so that the term Λ_γδ(${\boldsymbol{K}}$, ${\boldsymbol{K}}$') never contributes to Equation (66). For ℓ_α even, Equation (66) becomes

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}{C}_{\alpha \beta }({\boldsymbol{K}},{\boldsymbol{K}}^{\prime} ,t)\right|}_{t=0}=-{\delta }_{\alpha }^{\beta }\,{A}_{{{\ell }}_{\alpha }}\,{\delta }_{{\rm{D}}}({\boldsymbol{K}}-{\boldsymbol{K}}^{\prime} )\,n({\boldsymbol{K}})\,{{\rm{\Gamma }}}^{2}({\boldsymbol{K}}),\end{eqnarray} \tag{ 68 }$$

where the sum over ℓ_δ is performed following Equation (60), and the decay rate Γ² (${\boldsymbol{K}}$) is given by Equation (15).

In this appendix, we compute the two averages appearing in Equation (32). Assuming that the test particle is initially uniformly distributed on the sphere, one straightforwardly has

$$\begin{eqnarray}&&{\left\langle {\varphi }_{\alpha }^{{\rm{t}}}(0){\varphi }_{\beta }^{{\rm{t}}}(0)\right\rangle }_{{\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}}=\int \displaystyle \frac{d{\widehat{{\boldsymbol{L}}}}_{{\rm{t}}}}{4\pi }\,{Y}_{\alpha }({\widehat{{\boldsymbol{L}}}}_{{\rm{t}}})\,{Y}_{\beta }({\widehat{{\boldsymbol{L}}}}_{{\rm{t}}})=\displaystyle \frac{{\delta }_{\alpha }^{\beta }}{4\pi }.\end{eqnarray} \tag{ 69 }$$

To compute the time average of $e$^Ω (t), we rely on the cumulant theorem,

$$\begin{eqnarray}&&{\left\langle {e}^{\varepsilon {\rm{\Omega }}}\right\rangle }_{T}=\displaystyle \sum _{n=0}^{\infty }\displaystyle \frac{{\varepsilon }^{n}}{n!}\,{\mu }_{n}=\exp \left[\displaystyle \sum _{n=1}^{\infty }\displaystyle \frac{{\varepsilon }^{n}}{n!}{\kappa }_{n}\right],\end{eqnarray} \tag{ 70 }$$

where ${\mu }_{n}={\left\langle {{\rm{\Omega }}}^{n}\right\rangle }_{T}$ are the moment matrices, and κ_n are the cumulant matrices, with the first two given by κ₁ = μ₁ and ${\kappa }_{2}={\mu }_{2}-{\mu }_{1}^{2}$. Here, we compute the time average of $e$^Ω(t) by keeping only terms that are at most second order in ${Q}^{{\rm{t}}}$, so that only $\langle {{\rm{\Omega }}}_{1}(t){\rangle }_{T}\propto {Q}^{{\rm{t}}}$ and $\langle {{\rm{\Omega }}}_{1}^{2}{(t){\rangle }_{T}\propto \langle {{\rm{\Omega }}}_{2}(t){\rangle }_{T}\propto ({Q}^{{\rm{t}}})}^{2}$ contribute. We note that since $\langle {\varphi }_{\alpha }({\boldsymbol{K}},t){\rangle }_{T}=0$ and, from stationarity, $\langle {Q}^{{\rm{t}}}(t)\,{Q}^{{\rm{t}}}({t}^{{\prime} }){\rangle }_{T}=\langle {Q}^{{\rm{t}}}({t}^{{\prime} })\,{Q}^{{\rm{t}}}(t){\rangle }_{T}$, both $\langle {{\rm{\Omega }}}_{1}(t){\rangle }_{T}$ and $\langle {{\rm{\Omega }}}_{2}(t){\rangle }_{T}$ vanish. As a result, at the order considered here, only the second cumulant ${\kappa }_{2}=\langle {{\rm{\Omega }}}_{1}^{2}(t){\rangle }_{T}$ is non-zero, and from the cumulant theorem we obtain

$$\begin{eqnarray}&&{\left[{\left\langle {e}^{{\rm{\Omega }}(t)}\right\rangle }_{T}\right]}_{\alpha \beta }={\left[{e}^{\tfrac{1}{2}\langle {{\rm{\Omega }}}_{1}^{2}(t){\rangle }_{T}}\right]}_{\alpha \beta }.\end{eqnarray} \tag{ 71 }$$

Gathering Equations (69) and (71) together, we can write the test particle's correlation function as

$$\begin{eqnarray}&&{C}_{\alpha \beta }^{{\rm{t}},{\rm{r}}}(t)=\displaystyle \frac{1}{4\pi }\,{\left[{e}^{\tfrac{1}{2}\langle {{\rm{\Omega }}}_{1}^{2}(t){\rangle }_{T}}\right]}_{\alpha \beta }.\end{eqnarray} \tag{ 72 }$$

Using Equation (30), we can write

$$\begin{eqnarray}&&\begin{array}{l}{\left[{\left\langle {{\rm{\Omega }}}_{1}^{2}(t)\right\rangle }_{T}\right]}_{\alpha \beta }={\displaystyle \int }_{0}^{t}{{ds}}_{1}{\displaystyle \int }_{0}^{t}{{ds}}_{2}\displaystyle \int d{{\boldsymbol{K}}}_{1}\displaystyle \int d{{\boldsymbol{K}}}_{2}\,{E}_{\alpha \gamma \delta }{E}_{\delta \lambda \beta }\\ \ \ \times \ {{ \mathcal J }}_{{{\ell }}_{\gamma }}\left[{{\boldsymbol{K}}}_{{\rm{t}}},{{\boldsymbol{K}}}_{1}\right]\,{{ \mathcal J }}_{{{\ell }}_{\lambda }}\left[{{\boldsymbol{K}}}_{{\rm{t}}},{{\boldsymbol{K}}}_{2}\right]\,{\left\langle {\varphi }_{\gamma }({{\boldsymbol{K}}}_{1},{s}_{1}){\varphi }_{\lambda }({{\boldsymbol{K}}}_{2},{s}_{2})\right\rangle }_{T}.\end{array}\end{eqnarray} \tag{ 73 }$$

To pursue the calculation further, we may now use our ansatz for the time-dependence of the correlation of the noise fluctuations, as obtained in Equation (22). Using the sum identities from Equation (60), one gets

$$\begin{eqnarray}&&{\left[{\left\langle {{\rm{\Omega }}}_{1}^{2}(t)\right\rangle }_{T}\right]}_{\alpha \beta }=-{\delta }_{\alpha }^{\beta }\,{A}_{{{\ell }}_{\alpha }}\displaystyle \sum _{{\ell }}{B}_{{\ell }}\int d{\boldsymbol{K}}\,{{ \mathcal J }}_{{\ell }}^{2}\left[{{\boldsymbol{K}}}_{{\rm{t}}},{\boldsymbol{K}}\right]{\chi }_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}},t),\end{eqnarray} \tag{ 74 }$$

where ${\chi }_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}},t)$ is a double time integral of the noise correlation

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}},t) & = & {\displaystyle \int }_{0}^{t}{{dt}}_{1}{\displaystyle \int }_{0}^{t}{{dt}}_{2}\,{C}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}},{t}_{1}-{t}_{2}).\\ & = & \ {n}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}})\,\displaystyle \frac{2{T}_{{\rm{c}}}^{2}({\boldsymbol{K}})}{{A}_{{\ell }}}\chi \left[\displaystyle \frac{t\sqrt{{A}_{{\ell }}/2}}{{T}_{{\rm{c}}}({\boldsymbol{K}})}\right],\end{array}\end{eqnarray} \tag{ 75 }$$

with the dimensionless function χ (τ) defined in Equation (34). Since the matrix $\langle {{\rm{\Omega }}}_{1}^{2}(t){\rangle }_{T}$ is diagonal, one can straightforwardly compute its exponential, as required by Equation (71). This allows us to finally recast the correlation of the test particle's random motion from Equation (72) under the form of Equation (33).

In this appendix, we compute the ensemble-averaged variance of the amplitude of the density fluctuations, ${C}_{\alpha \beta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },0)$, as introduced in Equation (20). Our goal is to compute an expression of the form

$$\begin{eqnarray}&&\begin{array}{l}\left\langle {C}_{\alpha \beta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },0)\,{C}_{\gamma \delta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\gamma },{{\boldsymbol{K}}}_{\delta },0)\right\rangle =\displaystyle \frac{1}{{T}^{2}}{\displaystyle \int }_{0}^{T}{dt}\\ \ \ \ \times \ {\displaystyle \int }_{0}^{T}{{dt}}^{{\prime} }\,\left\langle {\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha },t)\,{\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta },t)\,{\varphi }_{\gamma }({{\boldsymbol{K}}}_{\gamma },{t}^{{\prime} })\,{\varphi }_{\delta }({{\boldsymbol{K}}}_{\delta },{t}^{{\prime} })\right\rangle .\end{array}\end{eqnarray} \tag{ 76 }$$

Because only even harmonics contribute to the interactions, we can limit ourselves to 2 ≤ ℓ_α, ℓ_γ even. As estimated in Equation (23), we know that for $| t-{t}^{{\prime} }| \gg {T}_{{\rm{c}}}$, the location of the background particles at time ${t}^{{\prime} }$ can be considered to be decorrelated from their locations at time t, up to the requirement of satisfying the system's global conservation constraints. It is fundamental to account for these global conservation constraints, as they introduce non-ergodic effects, preventing us from interchanging time- and ensemble-averages. Provided that these constraints are satisfied, in the two-dimensional time integral from Equation (76), one can note that the particles are uncorrelated between t and ${t}^{{\prime} }$ on a surface of size (T − ${T}_{{\rm{c}}}$)², while they are correlated on a surface of size T ${T}_{{\rm{c}}}$. As a result, as long as T ≫ ${T}_{{\rm{c}}}$ and as long as the conservation constraints are satisfied, particles can be considered as uncorrelated between time t and ${t}^{{\prime} }$, and therefore distributed uniformly over the sphere at these two times.

Let us now detail how one may carry out the average from Equation (76), in the presence of these constraints. At time t, the state of the system is fully characterized by the set of all fields ${\boldsymbol{\varphi }}\,=\,\{{\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha },t)\}$ and, similarly, at time ${t}^{{\prime} }$, the state of the system is fully characterized by ${{\boldsymbol{\varphi }}}^{{\prime} }$ = {${\varphi }_{\alpha }^{{\prime} }({{\boldsymbol{K}}}_{\alpha },t)$}. We may then use ${\boldsymbol{\varphi }}$ and ${{\boldsymbol{\varphi }}}^{{\prime} }$ as the random variables over which averages are carried out. Following Equation (12), we have

$$\begin{eqnarray}\begin{array}{rcl}\left\langle {\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha },0)\right\rangle & = & 0\,(\mathrm{for}\ {{\ell }}_{\alpha }\ne 0);\\ \left\langle {\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha },0)\,{\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta },0)\right\rangle & = & {\delta }_{\alpha }^{\beta }\,{\delta }_{{\rm{D}}}({{\boldsymbol{K}}}_{\alpha }-{{\boldsymbol{K}}}_{\beta })\,n({{\boldsymbol{K}}}_{\alpha }),\end{array}\end{eqnarray} \tag{ 77 }$$

and similarly for ${{\boldsymbol{\varphi }}}^{{\prime} }$. Placing ourselves within the Gaussian limit, we may then treat ${\boldsymbol{\varphi }}$ (resp. ${{\boldsymbol{\varphi }}}^{{\prime} }$) as uncorrelated Gaussian random fields that follow a Gaussian PDF F (${\boldsymbol{\varphi }}$) (resp. F (${{\boldsymbol{\varphi }}}^{{\prime} }$)), with a covariance following from Equation (77).

As emphasized above, the two fields ${\boldsymbol{\varphi }}$ and ${{\boldsymbol{\varphi }}}^{{\prime} }$ remain correlated one with another through global constraints. To shorten the notation, let us temporarily denote these constraints as ${\boldsymbol{\theta }}$ = ${\boldsymbol{\theta }}$(${\boldsymbol{\varphi }}$). In Equation (76), the average must then be carried out according to the joint PDF $F({\boldsymbol{\varphi }},{{\boldsymbol{\varphi }}}^{{\prime} })\,=F({\boldsymbol{\varphi }})\,F({{\boldsymbol{\varphi }}}^{{\prime} }\,| \,{\boldsymbol{\theta }}({\boldsymbol{\varphi }}))$. The conditional PDF of ${{\boldsymbol{\varphi }}}^{{\prime} }$ given the constraint ${\boldsymbol{\theta }}$(${\boldsymbol{\varphi }}$) follows from Bayes' theorem, and reads

$$\begin{eqnarray}&&\begin{array}{r}F({{\boldsymbol{\varphi }}}^{{\prime} }\,| \,{\boldsymbol{\theta }}({\boldsymbol{\varphi }}))=\displaystyle \frac{F({{\boldsymbol{\varphi }}}^{{\prime} },{\boldsymbol{\theta }}({\boldsymbol{\varphi }}))}{{F}_{{\boldsymbol{\theta }}}({\boldsymbol{\theta }}({\boldsymbol{\varphi }}))}\\ \,\,\,\,=\displaystyle \frac{F({{\boldsymbol{\varphi }}}^{{\prime} })\,{\delta }_{{\rm{D}}}({\boldsymbol{\theta }}({\boldsymbol{\varphi }})-{\boldsymbol{\theta }}({{\boldsymbol{\varphi }}}^{{\prime} }))}{{F}_{{\boldsymbol{\theta }}}({\boldsymbol{\theta }}({\boldsymbol{\varphi }}))},\end{array}\end{eqnarray} \tag{ 78 }$$

with ${F}_{{\boldsymbol{\theta }}}({\boldsymbol{\theta }})$ the PDF of the constraints ${\boldsymbol{\theta }}$. Therefore, we can write

$$\begin{eqnarray}&&F({\boldsymbol{\varphi }},{{\boldsymbol{\varphi }}}^{{\prime} })=\displaystyle \frac{F({\boldsymbol{\varphi }})\,F({{\boldsymbol{\varphi }}}^{{\prime} })\,{\delta }_{{\rm{D}}}({\boldsymbol{\theta }}({\boldsymbol{\varphi }})-{\boldsymbol{\theta }}({{\boldsymbol{\varphi }}}^{{\prime} }))}{{F}_{{\boldsymbol{\theta }}}({\boldsymbol{\theta }}({\boldsymbol{\varphi }}))}.\end{eqnarray} \tag{ 79 }$$

In that view, Equation (76) can be recast as

$$\begin{eqnarray}&&\begin{array}{l}\left\langle {C}_{\alpha \beta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },0)\,{C}_{\gamma \delta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\gamma },{{\boldsymbol{K}}}_{\delta },0)\right\rangle \\ \ =\ \displaystyle \int d{\boldsymbol{\theta }}d{{\boldsymbol{\theta }}}^{{\prime} }\,\displaystyle \frac{{\delta }_{{\rm{D}}}({\boldsymbol{\theta }}-{{\boldsymbol{\theta }}}^{{\prime} })}{{F}_{{\boldsymbol{\theta }}}({\boldsymbol{\theta }})}\,{{\mathbb{F}}}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },{\boldsymbol{\theta }})\,{{\mathbb{F}}}_{\gamma \delta }({{\boldsymbol{K}}}_{\gamma },{{\boldsymbol{K}}}_{\delta },{{\boldsymbol{\theta }}}^{{\prime} }),\end{array}\end{eqnarray} \tag{ 80 }$$

where we introduce

$$\begin{eqnarray}&&\begin{array}{l}{{\mathbb{F}}}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },{\boldsymbol{\theta }})\\ \ \equiv \ \displaystyle \int d{\boldsymbol{\varphi }}\,F({\boldsymbol{\varphi }})\,{\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha })\,{\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta })\,{\delta }_{{\rm{D}}}({\boldsymbol{\theta }}({\boldsymbol{\varphi }})-{\boldsymbol{\theta }})\\ \ =\ \left\langle {\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha })\,{\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta })\,{\delta }_{{\rm{D}}}({\boldsymbol{\theta }}({\boldsymbol{\varphi }})-{\boldsymbol{\theta }})\right\rangle ,\end{array}\end{eqnarray} \tag{ 81 }$$

with $\langle \cdot \rangle$ standing for the ensemble average where the fields $\varphi$ are drawn according to the Gaussian statistics of F(${\boldsymbol{\varphi }}$). Conveniently, in that form, Equation (80) allows us to carry out independently the averages over ${\boldsymbol{\varphi }}$ and ${{\boldsymbol{\varphi }}}^{{\prime} }$.

To proceed further, let us now detail the global conservation constraints that have to be satisfied throughout the system's evolution. There are three such constraints, namely the conservation of each particle's individual parameters (${\boldsymbol{\theta }}$₀), the conservation of the system's total angular momentum (${\boldsymbol{\theta }}$₁), and the conservation of the system's total energy (${\boldsymbol{\theta }}$₂). Luckily, these can all be expressed as simple functions of the fields ${\boldsymbol{\varphi }}$. They read

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{l}{\boldsymbol{K}}{\textstyle \text{-}}{\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{e}}{\rm{r}}{\rm{v}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}\,({{\boldsymbol{\theta }}}_{0}):\,{\varphi }_{(0,0)}({\boldsymbol{K}},t)={\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{t}}.,\\ {\rm{Angular}}\,{\rm{m}}{\rm{o}}{\rm{m}}{\rm{e}}{\rm{n}}{\rm{t}}{\rm{u}}{\rm{m}}\,{\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{e}}{\rm{r}}{\rm{v}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}\,({{\boldsymbol{\theta }}}_{1}):\\ \displaystyle \frac{1}{N}\displaystyle \int d{\boldsymbol{K}}\,L\left[{\boldsymbol{K}}\right]\,{\varphi }_{(1,m)}({\boldsymbol{K}},t)={\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{t}}.\,{\rm{f}}{\rm{o}}{\rm{r}}\,-1\leqslant m\leqslant 1,\\ {\rm{E}}{\rm{n}}{\rm{e}}{\rm{r}}{\rm{g}}{\rm{y}}\,{\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{e}}{\rm{r}}{\rm{v}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}\,({{\boldsymbol{\theta }}}_{2}):\\ E=\displaystyle \frac{1}{2N}\displaystyle \sum _{\alpha }\displaystyle \int d{\boldsymbol{K}}d{{\boldsymbol{K}}}^{{\rm{{\prime} }}}\,{H}_{{\ell }_{\alpha }}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}^{{\rm{{\prime} }}}\right]\,{\varphi }_{\alpha }({\boldsymbol{K}},t)\,{\varphi }_{\alpha }({{\boldsymbol{K}}}^{{\rm{{\prime} }}},t)={\rm{c}}{\rm{o}}{\rm{n}}{\rm{s}}{\rm{t}}.,\end{array}\end{array}\end{eqnarray} \tag{ 82 }$$

with $L[{\boldsymbol{K}}]=m\sqrt{{{GM}}_{\bullet }a(1-{e}^{2})}$ the norm of the angular momentum and H_ℓ[${\boldsymbol{K}}$, ${{\boldsymbol{K}}}^{{\prime} }$] = L [${\boldsymbol{K}}$] ${{ \mathcal J }}_{{\ell }}$[${\boldsymbol{K}}$, ${{\boldsymbol{K}}}^{{\prime} }$]. We also note the prefactor 1/(2N) in the definition of the energy that was introduced for later convenience. At this stage, it is important to note that each of these constraints involve different harmonics of the Gaussian random fields, namely ℓ = 0 for the conservation of ${\boldsymbol{K}}$, ℓ = 1 for the angular momentum, and 2 ≤ ℓ even for the energy. In the limit of Gaussian random fields, this implies that only the energy constraint contributes to a non-zero variance in Equation (80), as we will now argue.

Since only even harmonics contribute to the interactions, we can restrict ourselves to 2 ≤ ℓ_α even when computing ${{\mathbb{F}}}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },{\boldsymbol{\theta }})$. If ℓ_β = 0, 1, Equation (81) can be rewritten as

$$\begin{eqnarray}&&\begin{array}{l}{{\mathbb{F}}}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },{\boldsymbol{\theta }})\\ \ =\ \left\langle {\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta })\,{\delta }_{D}({{\boldsymbol{\theta }}}_{0}({\boldsymbol{\varphi }})-{{\boldsymbol{\theta }}}_{0}){\delta }_{{\rm{D}}}({{\boldsymbol{\theta }}}_{1}({\boldsymbol{\varphi }})-{{\boldsymbol{\theta }}}_{1})\right\rangle \\ \ \ \ \times \ \left\langle {\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha })\,{\delta }_{{\rm{D}}}({{\boldsymbol{\theta }}}_{2}({\boldsymbol{\varphi }})-{{\boldsymbol{\theta }}}_{2})\right\rangle ,\end{array}\end{eqnarray} \tag{ 83 }$$

where we used the Gaussian assumption, so that fields with different harmonics are uncorrelated. Because the energy is quadratic in the fields, and because the Gaussian PDF, F(${\boldsymbol{\varphi }}$_ℓ≥2), is an even function of the fields, the last bracket in Equation (83) is equal to zero. As a result, we can assume that ℓ_β ≥ 2. In that case, Equation (81) becomes

$$\begin{eqnarray}&&\begin{array}{l}{{\mathbb{F}}}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },{\boldsymbol{\theta }})\\ \ =\ \left\langle {\delta }_{{\rm{D}}}({{\boldsymbol{\theta }}}_{0}({\boldsymbol{\varphi }})-{{\boldsymbol{\theta }}}_{0})\right\rangle \,\left\langle {\delta }_{{\rm{D}}}({{\boldsymbol{\theta }}}_{1}({\boldsymbol{\varphi }})-{{\boldsymbol{\theta }}}_{1})\right\rangle \\ \ \ \ \times \ \left\langle {\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha })\,{\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta })\,{\delta }_{{\rm{D}}}({{\boldsymbol{\theta }}}_{2}({\boldsymbol{\varphi }})-{{\boldsymbol{\theta }}}_{2})\right\rangle \\ \ =\ {F}_{{{\boldsymbol{\theta }}}_{0}}({{\boldsymbol{\theta }}}_{0})\,{F}_{{{\boldsymbol{\theta }}}_{1}}({{\boldsymbol{\theta }}}_{1})\,{{\mathbb{F}}}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },E),\end{array}\end{eqnarray} \tag{ 84 }$$

where we introduce

$$\begin{eqnarray}&&{{\mathbb{F}}}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },E)\equiv \left\langle {\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha })\,{\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta })\,{\delta }_{{\rm{D}}}(E({\boldsymbol{\varphi }})-E)\right\rangle .\end{eqnarray} \tag{ 85 }$$

As a result, for ℓ_α, ℓ_β, ℓ_γ, ℓ_δ ≥ 2, this allows us to rewrite the required correlation from Equation (80) as

$$\begin{eqnarray}&&\begin{array}{l}\left\langle {C}_{\alpha \beta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },0)\,{C}_{\gamma \delta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\gamma },{{\boldsymbol{K}}}_{\delta },0)\right\rangle =\displaystyle \int {{dEdE}}^{{\prime} }\,\displaystyle \frac{{\delta }_{{\rm{D}}}(E-{E}^{{\prime} })}{{F}_{E}(E)}\\ \quad \ \times \ {{\mathbb{F}}}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },E)\,{{\mathbb{F}}}_{\gamma \delta }({{\boldsymbol{K}}}_{\gamma },{{\boldsymbol{K}}}_{\delta },{E}^{{\prime} }),\end{array}\end{eqnarray} \tag{ 86 }$$

where we get rid of all occurrences of the constraints ${\boldsymbol{\theta }}$₀ and ${\boldsymbol{\theta }}$₁, using the fact that their PDFs satisfy $\int d{{\boldsymbol{\theta }}}_{0}\,{F}_{{{\boldsymbol{\theta }}}_{0}}=\int d{{\boldsymbol{\theta }}}_{1}{F}_{{{\boldsymbol{\theta }}}_{1}}=1$. As a conclusion, in the limit of Gaussian random fields, only the constraint of total energy conservation contributes to the non-ergodic properties of the system. This is an important result of this calculation.

Let us now compute the ensemble average appearing in Equation (85). In the present Gaussian limit, we can rely on Novikov's theorem (Novikov 1965) to compute it.⁷ One obtains

$$\begin{eqnarray}&&\begin{array}{l}{{\mathbb{F}}}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },E)=\displaystyle \sum _{\mu }\displaystyle \int d{{\boldsymbol{K}}}_{\mu }\,\left\langle {\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha })\,{\varphi }_{\mu }({{\boldsymbol{K}}}_{\mu })\right\rangle \\ \qquad \times \ \left\langle \displaystyle \frac{\delta }{\delta {\varphi }_{\mu }({{\boldsymbol{K}}}_{\mu })}\left[{\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta })\,{\delta }_{{\rm{D}}}(E({\boldsymbol{\varphi }})-E)\right]\right\rangle ,\end{array}\end{eqnarray} \tag{ 87 }$$

where the first cumulant is absent because ℓ_α ≥ 2, so that $\langle {\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha })\rangle =0$, and only the second cumulant remains as the fields are assumed to be Gaussian. In Equation (87), the sum (resp. integral) over μ (resp. $d{{\boldsymbol{K}}}_{\mu }$) runs over all the fields. The functional gradient appearing in the last term can be computed as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\delta }{\delta {\varphi }_{\mu }({{\boldsymbol{K}}}_{\mu })}\left[{\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta })\,{\delta }_{{\rm{D}}}(E({\boldsymbol{\varphi }})-E)\right]\\ \ =\ {\delta }_{\mu }^{\beta }\,{\delta }_{{\rm{D}}}({{\boldsymbol{K}}}_{\mu }-{{\boldsymbol{K}}}_{\beta })\,{\delta }_{{\rm{D}}}(E({\boldsymbol{\varphi }})-E)\\ \ \ -\ {\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta })\,\displaystyle \frac{\partial }{\partial E}\left[\displaystyle \frac{\delta E({\boldsymbol{\varphi }})}{\delta {\varphi }_{\mu }({{\boldsymbol{K}}}_{\mu })}\,{\delta }_{{\rm{D}}}(E({\boldsymbol{\varphi }})-E)\right],\end{array}\end{eqnarray} \tag{ 88 }$$

where we use the fundamental relation $\delta {\varphi }_{\beta }$ $({{\boldsymbol{K}}}_{\beta })$/$\delta {\varphi }_{\mu }$ $({{\boldsymbol{K}}}_{\mu })$ = ${\delta }_{\mu }^{\beta }{\delta }_{{\rm{D}}}({{\boldsymbol{K}}}_{\beta }$ − K_μ). Glancing back at the definition of the energy in Equation (82), we can also write

$$\begin{eqnarray}&&\displaystyle \frac{\delta E({\boldsymbol{\varphi }})}{\delta {\varphi }_{\mu }({{\boldsymbol{K}}}_{\mu })}=\displaystyle \frac{1}{N}\int d{\boldsymbol{K}}\,{H}_{{{\ell }}_{\mu }}\left[{{\boldsymbol{K}}}_{\mu },{\boldsymbol{K}}\right]\,{\varphi }_{\mu }({\boldsymbol{K}}).\end{eqnarray} \tag{ 89 }$$

Injecting these results into Equation (87) and using the Gaussian statistics from Equation (77), we obtain a self-consistent integro-differential equation for ${\mathbb{F}}$_αβ$({{\boldsymbol{K}}}_{\alpha }$, K_β, E), namely

$$\begin{eqnarray}&&\begin{array}{l}{{\mathbb{F}}}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },E)={\delta }_{\alpha }^{\beta }\,{\delta }_{{\rm{D}}}({{\boldsymbol{K}}}_{\alpha }-{{\boldsymbol{K}}}_{\beta })\,n({{\boldsymbol{K}}}_{\alpha })\,{F}_{E}(E)\\ \ \ -\ \displaystyle \frac{1}{N}\,n({{\boldsymbol{K}}}_{\alpha })\displaystyle \int d{\boldsymbol{K}}\,{H}_{{{\ell }}_{\alpha }}\left[{{\boldsymbol{K}}}_{\alpha },{\boldsymbol{K}}\right]\,\displaystyle \frac{\partial }{\partial E}{{\mathbb{F}}}_{\alpha \beta }({\boldsymbol{K}},{{\boldsymbol{K}}}_{\beta },E),\end{array}\end{eqnarray} \tag{ 90 }$$

where we use that $\langle {\delta }_{{\rm{D}}}(E({\boldsymbol{\varphi }})-E)\rangle$ = F_E (E), by definition.

Progress can now be made by accounting perturbatively for the total energy constraint. As such, we introduce the small parameter $\varepsilon$, make the substitution H_ℓ → $\varepsilon$H_ℓ in Equation (90), and consider the expansion

$$\begin{eqnarray}&&{\mathbb{F}}={{\mathbb{F}}}^{(0)}+\varepsilon \,{{\mathbb{F}}}^{(1)}+{\varepsilon }^{2}\,{{\mathbb{F}}}^{(2)}+....\end{eqnarray} \tag{ 91 }$$

We can then inject this expansion in Equation (90) and match the orders in $\varepsilon$. The first three terms are obtained as

$$\begin{eqnarray}\begin{array}{rcl}{{\mathbb{F}}}_{\alpha \beta }^{(0)}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },E) & = & {\delta }_{\alpha }^{\beta }\,{\delta }_{{\rm{D}}}({{\boldsymbol{K}}}_{\alpha }-{{\boldsymbol{K}}}_{\beta })\,n({{\boldsymbol{K}}}_{\alpha })\,{F}_{E}(E),\\ {{\mathbb{F}}}_{\alpha \beta }^{(1)}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },E) & = & -\displaystyle \frac{1}{N}{\delta }_{\alpha }^{\beta }\,n({{\boldsymbol{K}}}_{\alpha })\,n({{\boldsymbol{K}}}_{\beta })\\ & & \times \ \displaystyle \frac{\partial {F}_{E}(E)}{\partial E}\,{H}_{{{\ell }}_{\alpha }}\left[{{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta }\right],\\ {{\mathbb{F}}}_{\alpha \beta }^{(2)}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },E) & = & \displaystyle \frac{1}{{N}^{2}}{\delta }_{\alpha }^{\beta }\,n({{\boldsymbol{K}}}_{\alpha })\,n({{\boldsymbol{K}}}_{\beta })\,\displaystyle \frac{{\partial }^{2}{F}_{E}(E)}{\partial {E}^{2}}\\ & & \times \ \displaystyle \int d{\boldsymbol{K}}\,{H}_{{{\ell }}_{\alpha }}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}_{\beta }\right]\,{H}_{{{\ell }}_{\alpha }}\left[{{\boldsymbol{K}}}_{\alpha },{\boldsymbol{K}}\right]\,n({\boldsymbol{K}}).\end{array}\end{eqnarray} \tag{ 92 }$$

Owing to these first terms, we can now return to the computation of the variance from Equation (86). Keeping only terms of at most second order in $\varepsilon$, this reads

$$\begin{eqnarray}&&\begin{array}{l}\langle {C}_{\alpha \beta }^{{\rm{r}}}\,{C}_{\gamma \delta }^{{\rm{r}}}\rangle \\ \ =\ \langle {C}_{\alpha \beta }^{{\rm{r}}}\,{C}_{\gamma \delta }^{{\rm{r}}}{\rangle }^{(0)}+\varepsilon \,\langle {C}_{\alpha \beta }^{{\rm{r}}}\,{C}_{\gamma \delta }^{{\rm{r}}}{\rangle }^{(1)}+{\varepsilon }^{2}\,\langle {C}_{\alpha \beta }^{{\rm{r}}}\,{C}_{\gamma \delta }^{{\rm{r}}}{\rangle }^{(2)}+...\\ \ =\ \displaystyle \int {{dEdE}}^{{\prime} }\,\displaystyle \frac{{\delta }_{{\rm{D}}}(E-{E}^{{\prime} })}{{F}_{E}(E)}\left[{{\mathbb{F}}}_{\alpha \beta }^{(0)}(E)\,{{\mathbb{F}}}_{\gamma \delta }^{(0)}({E}^{{\prime} })\right.\\ \ \ +\ \varepsilon \left\{{{\mathbb{F}}}_{\alpha \beta }^{(0)}(E)\,{{\mathbb{F}}}_{\gamma \delta }^{(1)}({E}^{{\prime} })+{{\mathbb{F}}}_{\alpha \beta }^{(1)}(E)\,{{\mathbb{F}}}_{\gamma \delta }^{(0)}({E}^{{\prime} })\right\}\\ \ \ +\ {\varepsilon }^{2}\left\{{{\mathbb{F}}}_{\alpha \beta }^{(0)}(E)\,{{\mathbb{F}}}_{\gamma \delta }^{(2)}({E}^{{\prime} })+{{\mathbb{F}}}_{\alpha \beta }^{(1)}(E)\,{{\mathbb{F}}}_{\gamma \delta }^{(1)}({E}^{{\prime} })\right.\\ \ \ +\ \left.\left.{{\mathbb{F}}}_{\alpha \beta }^{(2)}(E)\,{{\mathbb{F}}}_{\gamma \delta }^{(0)}({E}^{{\prime} })\right\}\right],\end{array}\end{eqnarray} \tag{ 93 }$$

where, for simplicity, we do not repeat the arguments $({{\boldsymbol{K}}}_{\alpha }$, K_β) and $({{\boldsymbol{K}}}_{\gamma }$, K_δ). The zeroth-order term is straightforward to compute, and gives

$$\begin{eqnarray}&&\begin{array}{l}\langle {C}_{\alpha \beta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },0)\,{C}_{\gamma \delta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\gamma },{{\boldsymbol{K}}}_{\delta },0){\rangle }^{(0)}\\ \ =\ {\delta }_{\alpha }^{\beta }\,{\delta }_{\gamma }^{\delta }\,{\delta }_{{\rm{D}}}({{\boldsymbol{K}}}_{\alpha }-{{\boldsymbol{K}}}_{\beta })\,{\delta }_{{\rm{D}}}({{\boldsymbol{K}}}_{\gamma }-{{\boldsymbol{K}}}_{\delta })\,n({{\boldsymbol{K}}}_{\alpha })\,n({{\boldsymbol{K}}}_{\gamma })\\ \ =\ {C}_{\alpha \beta }({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },0)\,{C}_{\gamma \delta }({{\boldsymbol{K}}}_{\gamma },{{\boldsymbol{K}}}_{\delta },0),\end{array}\end{eqnarray} \tag{ 94 }$$

using $\int {dE}\,{F}_{E}(E)=1$. It is straighforward to show that the first-order term satisfies

$$\begin{eqnarray}&&\langle {C}_{\alpha \beta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },0)\,{C}_{\gamma \delta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\gamma },{{\boldsymbol{K}}}_{\delta },0){\rangle }^{(1)}=0,\end{eqnarray} \tag{ 95 }$$

as the energy integrals vanish. Using a similar argument, one finds that terms of the form ${\mathbb{F}}$⁽⁰⁾ ${\mathbb{F}}$⁽²⁾ do not contribute to the second-order term in Equation (93). Keeping only the non-zero contribution coming from ${{\mathbb{F}}}_{\alpha \beta }^{(1)}(E)\,{{\mathbb{F}}}_{\gamma \delta }^{(1)}({E}^{{\prime} })$, we get

$$\begin{eqnarray}&&\begin{array}{l}\langle {C}_{\alpha \beta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta },0)\,{C}_{\gamma \delta }^{{\rm{r}}}({{\boldsymbol{K}}}_{\gamma },{{\boldsymbol{K}}}_{\delta },0){\rangle }^{(2)}\\ \ =\ \displaystyle \frac{1}{{N}^{2}}\,{\delta }_{\alpha }^{\beta }\,{\delta }_{\gamma }^{\delta }\,n({{\boldsymbol{K}}}_{\alpha })\,n({{\boldsymbol{K}}}_{\beta })\,n({{\boldsymbol{K}}}_{\gamma })\,n({{\boldsymbol{K}}}_{\delta })\\ \ \ \times \ \displaystyle \frac{{H}_{{{\ell }}_{\alpha }}\left[{{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\beta }\right]\,{H}_{{{\ell }}_{\gamma }}\left[{{\boldsymbol{K}}}_{\gamma },{{\boldsymbol{K}}}_{\delta }\right]}{{\left({\rm{\Delta }}E\right)}^{2}},\end{array}\end{eqnarray} \tag{ 96 }$$

where ΔE is obtained after straightforward manipulation of the energy integrals, and reads

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\left({\rm{\Delta }}E\right)}^{2}}\equiv \int {dE}\,\displaystyle \frac{1}{{F}_{E}(E)}{\left(\displaystyle \frac{\partial {F}_{E}(E)}{\partial E}\right)}^{2}.\end{eqnarray} \tag{ 97 }$$

It is important to note that ΔE is a single number that depends only on the total energy PDF, F_E (E), and therefore only on the considered DF n(${\boldsymbol{K}}$). In Appendix G, we detail how the required integral from Equation (97) can be estimated. Equation (96) is an important result of this appendix, as it characterizes the variance (over different realizations) of the noise fluctuations' amplitude arising from the constraint of total energy conservation.

Following Equation (24), we can get the variance of ${n}_{{\ell }}^{r}$(${\boldsymbol{K}}$). It reads

$$\begin{eqnarray}&&\left\langle {n}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}})\,{n}_{{{\ell }}^{{\prime} }}^{{\rm{r}}}({{\boldsymbol{K}}}^{{\prime} })\right\rangle =n({\boldsymbol{K}})\,n({{\boldsymbol{K}}}^{{\prime} })\left\{1+{M}_{{\ell }}({\boldsymbol{K}})\,{M}_{{{\ell }}^{{\prime} }}({{\boldsymbol{K}}}^{{\prime} })\right\},\end{eqnarray} \tag{ 98 }$$

where we introduce the dimensionless function M_ℓ(${\boldsymbol{K}}$) as

$$\begin{eqnarray}&&{M}_{{\ell }}({\boldsymbol{K}})=\int d{{\boldsymbol{K}}}^{{\rm{{\prime} }}}\,\displaystyle \frac{n({{\boldsymbol{K}}}^{{\prime} })}{N}\,\displaystyle \frac{{H}_{{\ell }}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} }\right]}{{\rm{\Delta }}E}.\end{eqnarray} \tag{ 99 }$$

The final step of this appendix is to compute the variance of ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$, as defined in Equation (37). We get

$$\begin{eqnarray}&&\langle {\left({\rm{\Delta }}{{\rm{\Gamma }}}_{{\rm{t}}}^{2}\right)}^{2}\rangle ={\left\{\displaystyle \sum _{{\ell }}{B}_{{\ell }}\int d{\boldsymbol{K}}n({\boldsymbol{K}}){{ \mathcal J }}_{{\ell }}^{2}\left[{{\boldsymbol{K}}}_{{\rm{t}}},{\boldsymbol{K}}\right]{M}_{{\ell }}({\boldsymbol{K}})\right\}}^{2},\end{eqnarray} \tag{ 100 }$$

with ${\rm{\Delta }}{{\rm{\Gamma }}}_{{\rm{t}}}^{2}$ = ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$ − Γ². Equation (100) is the final result of this section. It expresses the variance (over realizations) of ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$, the amplitude of the random walk of a given test star. It is important to note that since (Γ²)² and $\langle {({\rm{\Delta }}{{\rm{\Gamma }}}_{{\rm{t}}}^{2})}^{2}\rangle$ have the same scaling with N, the present variance effect does not vanish as the number of particles gets larger.

As it will be needed to obtain the prediction of Figure 6, let us illustrate the effect associated with this non-zero variance of ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$ in our fiducial simulations. We consider the same window, ${W}_{{\rm{t}}}$(${{\boldsymbol{K}}}_{{\rm{t}}}$), as in Equation (115). For each test particle falling in that window, we measure the correlation function $\langle {\varphi }_{\alpha }^{{\rm{t}}}(t)\,{\varphi }_{\alpha }^{{\rm{t}}}(0){\rangle }_{T}$ (e.g., for ℓ_α = 1). The second-order time derivative at t = 0 of this correlation function is directly proportional to ${{\rm{\Gamma }}}_{{\rm{t}},\mathrm{num}.}^{2}$ (see Equation (36)), which we can therefore measure numerically. In Figure 8, we represent the distribution of these numerically measured initial values, ${{\rm{\Gamma }}}_{{\rm{t}},\mathrm{num}.}^{2}$, and illustrate how these amplitudes vary from realization to realization.

Zoom In Zoom Out Reset image size

To capture this variance effect seen in Figure 8, we may use our estimation of the variance of ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$ obtained in Equation (100). To do so, for every test particle falling in the window, we compute $({{\rm{\Gamma }}}^{2},\langle {({\rm{\Delta }}{{\rm{\Gamma }}}_{{\rm{t}}}^{2})}^{2}\rangle )$, following Equations (38) and (100). Having determined this mean and variance, one can draw a sample of ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$ according to a PDF that shares the same first two cumulants. Similarly to Equation (105), we assume that we can draw this sample according to a chi-squared PDF with these imposed mean and variance. The result of this procedure is illustrated in Figure 8. This figure illustrates how the present calculations are able to capture most of the features associated with the non-vanishing variance of the test particles' individual amplitudes ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$.

In this appendix, we briefly detail how one can estimate the energy spread ΔE introduced in Equation (97). Our convention for the definition of the total energy is spelled out in Equation (82). Using the two-point statistics from Equation (62), the expected mean value of the energy reads

$$\begin{eqnarray}&&\left\langle E\right\rangle =\displaystyle \frac{1}{2}\displaystyle \sum _{{\ell }}(2{\ell }+1)\,{E}_{{\ell }}^{(1)},\end{eqnarray} \tag{ 101 }$$

where we define

$$\begin{eqnarray}&&{E}_{{\ell }}^{(1)}=\int d{\boldsymbol{K}}\,{H}_{{\ell }}\left[{\boldsymbol{K}},{\boldsymbol{K}}\right]\,\displaystyle \frac{n({\boldsymbol{K}})}{N}.\end{eqnarray} \tag{ 102 }$$

We can proceed similarly to compute the expectation for $\langle {E}^{2}\rangle$. This reads

$$\begin{eqnarray}&&\begin{array}{l}\langle {E}^{2}\rangle =\displaystyle \frac{1}{4}\displaystyle \sum _{\alpha ,\beta }\displaystyle \int d{{\boldsymbol{K}}}_{\alpha }{{\boldsymbol{K}}}_{\alpha }^{{\prime} }d{{\boldsymbol{K}}}_{\beta }d{{\boldsymbol{K}}}_{\beta }^{{\prime} }\,{H}_{{{\ell }}_{\alpha }}\left[{{\boldsymbol{K}}}_{\alpha },{{\boldsymbol{K}}}_{\alpha }^{{\prime} }\right]\\ \,\,\times \ {H}_{{{\ell }}_{\beta }}\left[{{\boldsymbol{K}}}_{\beta },{{\boldsymbol{K}}}_{\beta }^{{\prime} }\right]\,\left\langle {\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha })\,{\varphi }_{\alpha }({{\boldsymbol{K}}}_{\alpha }^{{\prime} })\,{\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta })\,{\varphi }_{\beta }({{\boldsymbol{K}}}_{\beta }^{{\prime} })\right\rangle .\end{array}\end{eqnarray} \tag{ 103 }$$

To compute the average term appearing in the rhs, we follow Equation (63), placing ourselves in the limit of Gaussian random fields so that only connected averages involving two fields remain. We get

$$\begin{eqnarray}&&\langle {E}^{2}\rangle -\langle E{\rangle }^{2}=\displaystyle \frac{1}{2}\displaystyle \sum _{{\ell }}(2{\ell }+1)\,{E}_{{\ell }}^{(2)},\end{eqnarray} \tag{ 104 }$$

where we use the symmetry relation H_ℓ[${\boldsymbol{K}}$, ${{\boldsymbol{K}}}^{{\prime} }$] = H_ℓ[${{\boldsymbol{K}}}^{{\prime} }$, ${\boldsymbol{K}}$], and introduce

$$\begin{eqnarray}&&{E}_{{\ell }}^{(2)}=\int d{\boldsymbol{K}}d{{\boldsymbol{K}}}^{{\rm{{\prime} }}}\,{H}_{{\ell }}^{2}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} }\right]\,\displaystyle \frac{n({\boldsymbol{K}})}{N}\,\displaystyle \frac{n({{\boldsymbol{K}}}^{{\prime} })}{N}.\end{eqnarray} \tag{ 105 }$$

Having estimated the mean and the variance of the energy distribution, we may now return to the evaluation of the energy spread ΔE introduced in Equation (97). As the energy is a quadratic function of Gaussian fields, we will assume that it follows a (scaled) chi-squared distribution of mean μ = $\langle {\text{}}E\rangle$, and variance ${\sigma }^{2}=\langle {E}^{2}\rangle -\langle E{\rangle }^{2}$. The associated PDF then follows

$$\begin{eqnarray}&&{F}_{E}(E)=\displaystyle \frac{1}{\mu }\displaystyle \frac{\kappa }{{\rm{\Gamma }}(\kappa )}{\left(\kappa \displaystyle \frac{E}{\mu }\right)}^{\kappa -1}{e}^{-\kappa E/\mu }.\end{eqnarray} \tag{ 106 }$$

with κ = μ²/σ². For our fiducial numerical system, we find $\langle {\text{}}E\rangle$ ≃ 6.1 × 10⁻³ and κ ≃ 81. In Figure 9, we illustrate the statistical distribution of the system's energy, as well as the approximation from Equation (106).

Zoom In Zoom Out Reset image size

Finally, for a chi-squared PDF as in Equation (106), one can explicitly compute ΔE, as defined in Equation (97), to get

$$\begin{eqnarray}&&{\rm{\Delta }}E=\left\langle E\right\rangle \displaystyle \frac{\sqrt{\kappa -2}}{\kappa }.\end{eqnarray} \tag{ 107 }$$

We note that this integral is well-behaved only for κ > 2. This is an artefact coming from the perturbative expansion introduced in Equation (91). For our fiducial model, we find ΔE ≃ 6.7 × 10⁻⁴.

In this appendix, we briefly detail the procedures used in Figures 3 and 6 to compare our analytical results with the fiducial numerical simulations.

H.1. Correlation of the Noise Fluctuations

Let us detail the method followed to obtain Figure 3, used to illustrate Equation (22). One of the key insight from this equation is that to any particle (of parameter ${\boldsymbol{K}}$), we can associate the pair (n(${\boldsymbol{K}}$), ${T}_{{\rm{c}}}$(${\boldsymbol{K}}$)) that characterizes the correlation properties of the density fluctuations generated by background particles with these parameters. As a result, in order to consider only particles that have similar noise decorrelation properties, it is convenient to introduce, for every realization, the ${\boldsymbol{K}}$-averaged fields ${\overline{\varphi }}_{\alpha }(t)\,=\int d{\boldsymbol{K}}\,W({\boldsymbol{K}})\,{\varphi }_{\alpha }({\boldsymbol{K}},t)$, with W(${\boldsymbol{K}}$) a window function defined⁸ as

$$\begin{eqnarray}W({\boldsymbol{K}})=\left\{\begin{array}{cc}1 & \mathrm{if}\ 1\leqslant \displaystyle \frac{n({\boldsymbol{K}})}{{C}_{\min }},\displaystyle \frac{{T}_{{\rm{c}}}({\boldsymbol{K}})}{{T}_{{\rm{c}}}^{\min }}\leqslant 1+{\varepsilon }_{W},\\ 0 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 108 }$$

with (C_min, ${T}_{{\rm{c}}}^{\min }$) the typical amplitude and torque time considered, and $\varepsilon$_W a small dimensionless parameter controlling the size of the window. We then naturally have ${\overline{\varphi }}_{\alpha }(t)\,={\sum }_{i\in W}{Y}_{\alpha }({\widehat{{\boldsymbol{L}}}}_{i}(t))$, with the sum limited to the particles with (n(${\boldsymbol{K}}$), ${T}_{{\rm{c}}}$(${\boldsymbol{K}}$)) in the vicinity of (C_min, ${T}_{{\rm{c}}}^{\min }$). Introducing ${C}_{{{\ell }}_{\alpha },W}^{{\rm{r}}}$(t − ${t}^{{\prime} }$) ≡ $\langle {\overline{\varphi }}_{\alpha }(t)\,{\overline{\varphi }}_{\alpha }({t}^{{\prime} }){\rangle }_{T}$, Equation (22) immediately gives

$$\begin{eqnarray}&&{C}_{{\ell },W}^{{\rm{r}}}(t)=\int d{\boldsymbol{K}}\,W({\boldsymbol{K}})\,{C}_{{\ell }}^{{\rm{r}}}({\boldsymbol{K}},t),\end{eqnarray} \tag{ 109 }$$

where the function ${C}_{{\ell }}^{{\rm{r}}}$(${\boldsymbol{K}}$, t) follows the Gaussian ansatz from Equation (23). When averaged over realizations, Equation (109) can be approximated with the Gaussian dependence

$$\begin{eqnarray}&&{C}_{{\ell },W}^{{\rm{r}}}(t)\simeq {C}_{W}\,{e}^{-\tfrac{{A}_{{\ell }}}{2}{\left(t/{T}_{W}\right)}^{2}},\end{eqnarray} \tag{ 110 }$$

where we introduce the amplitude C_W and torque time T_W as

$$\begin{eqnarray}&&{C}_{W}=\int d{\boldsymbol{K}}\,W({\boldsymbol{K}})\,n({\boldsymbol{K}});\,\displaystyle \frac{{C}_{W}}{{\left({T}_{W}\right)}^{2}}=\int d{\boldsymbol{K}}\,W({\boldsymbol{K}})\displaystyle \frac{n({\boldsymbol{K}})}{{\left({T}_{{\rm{c}}}({\boldsymbol{K}})\right)}^{2}}.\end{eqnarray} \tag{ 111 }$$

Equation (110) is the analytical Gaussian prediction illustrated in Figure 3.

In Figure 3, we also present an updated prediction of the noise correlation obtained by reinjecting the Gaussian prediction from Equation (18) into the self-consistency relation from Equation (44). In that context, when averaged over the window, the prediction takes the form

$$\begin{eqnarray}&&{C}_{\ell ,W}(t)\simeq {C}_{W}\,\,\exp \left\{-\displaystyle \frac{{A}_{\ell }}{2}{\int }_{0}^{t}{dt}_{1}{\int }_{0}^{t}{dt}_{2}\,{{\rm{\Gamma }}}_{W}^{2}\,{\exp }^{-{\left[\left({t}_{1}-{t}_{2})/{T}_{W}^{{\rm{t}}}\right.\right]}^{2}}\right\},\end{eqnarray} \tag{ 112 }$$

with C_W as in Equation (111), and where the amplitude ${{\rm{\Gamma }}}_{W}^{2}$ and coherence time ${T}_{W}^{{\rm{t}}}$ are given by

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{W}^{2}={\left\langle {{\rm{\Gamma }}}_{{\rm{t}}}^{2}({\boldsymbol{K}})\right\rangle }_{W};\,{T}_{W}^{{\rm{t}}}={\left\langle {T}_{{\rm{c}}}^{{\rm{t}}}({\boldsymbol{K}})\right\rangle }_{W}.\end{eqnarray} \tag{ 113 }$$

In these equations, we introduce ${\left\langle \cdot \right\rangle }_{W}$ as the mean over the window W(${\boldsymbol{K}}$), i.e., it is defined as

$$\begin{eqnarray}&&{\left\langle f({\boldsymbol{K}})\right\rangle }_{W}=\displaystyle \frac{\int d{\boldsymbol{K}}\,n({\boldsymbol{K}})\,W({\boldsymbol{K}})\,f({\boldsymbol{K}})}{\int d{\boldsymbol{K}}\,n({\boldsymbol{K}})\,W({\boldsymbol{K}})}.\end{eqnarray} \tag{ 114 }$$

We illustrate such measurements in Figure 3. As defined in Equation (108), we chose the window W(${\boldsymbol{K}}$) with the parameters (C_min, ${T}_{{\rm{c}}}^{\min }$) ≃ (5.0, 142), with $\varepsilon$_W = 0.1. With such a choice, there are on average nine particles in the window per simulation. Following Equation (111), the typical amplitude and torque time in the window were found to be (C_W, T_W) ≃ (0.72, 148).

H.2. Correlation of the Random Walks

Let us briefly detail the method followed to obtain Figure 6, used to illustrate the result from Equation (36). Following the independence hypothesis from Equation (32) we assume that, for a given realization, each individual particle can effectively be treated as a test particle. As such, we neglect the correlations existing between the background fluctuations and the random walk of that particular particle.

One important insight from Equation (36) is that to any test particle (of parameter ${{\boldsymbol{K}}}_{{\rm{t}}}$), we can associate the pair (Γ²(${{\boldsymbol{K}}}_{{\rm{t}}}$), ${T}_{{\rm{c}}}^{{\rm{t}}}({{\boldsymbol{K}}}_{t}))$, that characterizes the correlation properties of its random walk in orientation. Similarly to Equation (108), in order to investigate these random walks, it is convenient to introduce the window function ${W}_{{\rm{t}}}$(${{\boldsymbol{K}}}_{{\rm{t}}}$) as

$$\begin{eqnarray}{W}_{{\rm{t}}}({{\boldsymbol{K}}}_{{\rm{t}}})=\left\{\begin{array}{cc}1 & \mathrm{if}\ 1\leqslant \displaystyle \frac{{{\rm{\Gamma }}}^{2}({{\boldsymbol{K}}}_{{\rm{t}}})}{{{\rm{\Gamma }}}_{\min }^{2}},\displaystyle \frac{{T}_{{\rm{c}}}^{{\rm{t}}}({{\boldsymbol{K}}}_{{\rm{t}}})}{{T}_{{\rm{t}}}^{\min }}\leqslant 1+{\varepsilon }_{{W}_{{\rm{t}}}},\\ 0 & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 115 }$$

where (${{\rm{\Gamma }}}_{\min }^{2}$, ${T}_{{\rm{t}}}^{\min }$) are the typical amplitude and coherence time of the considered test particles, and ${\varepsilon }_{{W}_{{\rm{t}}}}$ is a dimensionless parameter controlling the size of the window. We may then define the window-averaged correlation function (which can easily be measured in the numerical simulations) as

$$\begin{eqnarray}&&{C}_{{{\ell }}_{\alpha },{W}_{{\rm{t}}}}^{{\rm{t}}}(t)=4\pi \left\langle {\left\langle {\varphi }_{\alpha }^{{\rm{t}}}(t){\varphi }_{\alpha }^{{\rm{t}}}(0)\right\rangle }_{{W}_{{\rm{t}}}}\right\rangle ,\end{eqnarray} \tag{ 116 }$$

where $\langle \cdot {\rangle }_{{W}_{{\rm{t}}}}$ stands for the mean over the test particles of a given realization falling in the window ${W}_{{\rm{t}}}$(${{\boldsymbol{K}}}_{{\rm{t}}}$), similarly to Equation (114). Following Equation (36), we also add a prefactor 4π to ensure that this correlation is between 0 and 1.

If one does not account for the variance in ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$ (see Equation (100)), a first (naive) prediction for Equation (116) can be obtained from Equation (36) by restricting ourselves only to the ensemble-averaged mean predictions. This gives

$$\begin{eqnarray}&&{C}_{{\ell },{W}_{{\rm{t}}}}^{{\rm{t}}}(t)\simeq \exp \left\{-\displaystyle \frac{{A}_{{\ell }}}{2}\,{\int }_{0}^{t}{{dt}}_{1}{\int }_{0}^{t}{{dt}}_{2}\,{{\rm{\Gamma }}}_{{W}_{{\rm{t}}}}^{2}\,{e}^{-{\left[({t}_{1}-{t}_{2})/{T}_{{W}_{{\rm{t}}}}^{{\rm{t}}}\right]}^{2}}\right\},\end{eqnarray} \tag{ 117 }$$

where the amplitude, ${{\rm{\Gamma }}}_{{W}_{{\rm{t}}}}^{2}$, and coherence time, ${T}_{{W}_{{\rm{t}}}}$, are computed by direct averages over the particles falling in the window ${W}_{{\rm{t}}}$(${{\boldsymbol{K}}}_{{\rm{t}}}$), so that

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{W}_{{\rm{t}}}}^{2}={\left\langle {{\rm{\Gamma }}}^{2}({{\boldsymbol{K}}}_{{\rm{t}}})\right\rangle }_{{W}_{{\rm{t}}}};\,{T}_{{W}_{{\rm{t}}}}^{{\rm{t}}}={\left\langle {T}_{{\rm{c}}}^{{\rm{t}}}({{\boldsymbol{K}}}_{{\rm{t}}})\right\rangle }_{{W}_{{\rm{t}}}}.\end{eqnarray} \tag{ 118 }$$

One can improve the prediction from Equation (117) by accounting for the variance in ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$. To do so, for every test particle falling in the window, one can compute the mean expectation for the amplitude, $\langle {{\rm{\Gamma }}}_{{\rm{t}}}^{2}({{\boldsymbol{K}}}_{{\rm{t}}})\rangle ={{\rm{\Gamma }}}^{2}({{\boldsymbol{K}}}_{{\rm{t}}})$, and the associated variance, $\langle {({\rm{\Delta }}{{\rm{\Gamma }}}_{{\rm{t}}}^{2}({{\boldsymbol{K}}}_{{\rm{t}}}))}^{2}\rangle$, as given by Equations (38) and (100). For this same particle, one can then draw an effective value of ${{\rm{\Gamma }}}_{{\rm{t}}}^{2}$, according to a chosen PDF with this prescribed mean and variance. This process is illustrated in Figure 8, where we used a chi-squared PDF. In that case, the prediction from Equation (117) becomes

$$\begin{eqnarray}&&{C}_{{\ell },{W}_{{\rm{t}}}}^{{\rm{t}}}(t)\simeq {\left\langle \exp \left\{-\displaystyle \frac{{A}_{{\ell }}}{2}{\int }_{0}^{t}{{dt}}_{1}{\int }_{0}^{t}{{dt}}_{2}{{\rm{\Gamma }}}_{{\rm{t}}}^{2}{e}^{-{\left[({t}_{1}-{t}_{2})/{T}_{{W}_{{\rm{t}}}}^{{\rm{t}}}\right]}^{2}}\right\}\right\rangle }_{{W}_{{\rm{t}}}}.\end{eqnarray} \tag{ 119 }$$

Here, the amplitude and coherence time from Equation (118) become

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{\rm{t}}}^{2}=\mathrm{PDF}\left[{{\rm{\Gamma }}}^{2},\langle {\left({\rm{\Delta }}{{\rm{\Gamma }}}_{{\rm{t}}}^{2}\right)}^{2}\rangle \right];\,{T}_{{W}_{{\rm{t}}}}^{{\rm{t}}}={T}_{{\rm{c}}}^{{\rm{t}}},\end{eqnarray} \tag{ 120 }$$

where PDF [μ, σ²] returns a sample from a chosen PDF of mean μ and variance σ² (which we chose to be a chi-squared PDF as in Figure 8).

We illustrate such measurements in Figure 6. As defined by Equation (115), we chose the window ${W}_{{\rm{t}}}$(${{\boldsymbol{K}}}_{{\rm{t}}}$) with the parameters $({{\rm{\Gamma }}}_{\min }^{2},{T}_{{\rm{t}}}^{\min })=(1.2\times {10}^{-4},40.0)$ and ${\varepsilon }_{{W}_{{\rm{t}}}}=0.1$. Following Equation (118), the typical amplitude and coherence time in the window were found to be $({{\rm{\Gamma }}}_{{W}_{{\rm{t}}}}^{2},{T}_{{W}_{{\rm{t}}}}^{{\rm{t}}})\,\simeq (1.3\times {10}^{-4},43.6)$.

In this appendix, we detail all the calculations presented in Section 6 for an infinite power-law stellar distribution around an MBH. We first note that the squared coupling coefficients from Equation (51) can be rewritten as

$$\begin{eqnarray}&&{{ \mathcal J }}_{{\ell }}^{2}\left[{\boldsymbol{K}},{{\boldsymbol{K}}}^{{\prime} }\right]=\displaystyle \frac{G}{{M}_{\bullet }}\,\displaystyle \frac{1}{a(1-{e}^{2})}\,\displaystyle \frac{{m}^{{\prime} 2}}{{a}_{\mathrm{out}}^{2}}\,{s}_{{\ell }}^{2}\left[\alpha ,{e}_{\mathrm{in}},{e}_{\mathrm{out}}\right],\end{eqnarray} \tag{ 121 }$$

where we recall that "out" (resp. "in") labels the star with the larger (resp. smaller) semimajor axis, and we introduced the dimensionless ratio α = ${a}_{\mathrm{in}}$/${a}_{\mathrm{out}}$.

Following Equation (15), we can then compute the amplitude, Γ², of the background fluctuations. It reads

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}^{2}({\boldsymbol{K}}) & = & \displaystyle \frac{1}{4\pi }\displaystyle \frac{G\,\langle {m}^{2}\rangle }{{M}_{\bullet }}\,\displaystyle \frac{1}{a(1-{e}^{2})}\displaystyle \sum _{{\ell }}{B}_{{\ell }}\\ & & \times \ \displaystyle \int {{da}}^{{\prime} }{{de}}^{{\prime} }\,\displaystyle \frac{{f}_{e}({e}^{{\prime} })\,{n}_{a}({a}^{{\prime} })}{{a}_{\mathrm{out}}^{2}}\,{s}_{{\ell }}^{2}\left[\tfrac{{a}_{\mathrm{in}}}{{a}_{\mathrm{out}}},{e}_{\mathrm{in}},{e}_{\mathrm{out}}\right],\end{array}\end{eqnarray} \tag{ 122 }$$

where we introduce the second moment of the mass distribution $\langle {m}^{2} \rangle =\int {dmf}_{m}(m)\,{m}^{2}$. The integral over ${a}^{{\prime} }$ in Equation (122) can then be split into two regions, ${a}^{{\prime} }$ ≤ a and ${a}^{{\prime} }$ ≥ a. For the first region, we write

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{0}^{a}{{da}}^{{\prime} }\,\displaystyle \frac{{n}_{a}({a}^{{\prime} })}{{a}_{\mathrm{out}}^{2}}\,{s}_{{\ell }}^{2}\left[\tfrac{{a}_{\mathrm{in}}}{{a}_{\mathrm{out}}},{e}_{\mathrm{in}},{e}_{\mathrm{out}}\right]\\ \quad \ =\ {\displaystyle \int }_{0}^{a}{{da}}^{{\prime} }\,\displaystyle \frac{{n}_{a}({a}^{{\prime} })}{{a}^{2}}\,{s}_{{\ell }}^{2}\left[\tfrac{{a}^{{\prime} }}{a},{e}^{{\prime} },e\right]\\ \quad \ =\ {\displaystyle \int }_{0}^{1}d\alpha \,\displaystyle \frac{{n}_{a}(a\alpha )}{a}\,{s}_{{\ell }}^{2}\left[\alpha ,{e}^{{\prime} },e\right]\\ \quad \ =\ \displaystyle \frac{{N}_{0}}{{a}_{0}^{2}}{\left(\displaystyle \frac{a}{{a}_{0}}\right)}^{1-\gamma }{\displaystyle \int }_{0}^{1}d\alpha \,{\alpha }^{2-\gamma }\,{s}_{{\ell }}^{2}\left[\alpha ,{e}^{{\prime} },e\right],\end{array}\end{eqnarray} \tag{ 123 }$$

and a very similar calculation can be carried out for the second region to get

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{a}^{+\infty }{{da}}^{{\prime} }\,\displaystyle \frac{{n}_{a}({a}^{{\prime} })}{{a}_{\mathrm{out}}^{2}}\,{s}_{{\ell }}^{2}\left[\tfrac{{a}_{\mathrm{in}}}{{a}_{\mathrm{out}}},{e}_{\mathrm{in}},{e}_{\mathrm{out}}\right]\\ \quad \ =\ \displaystyle \frac{{N}_{0}}{{a}_{0}^{2}}{\left(\displaystyle \frac{a}{{a}_{0}}\right)}^{1-\gamma }{\displaystyle \int }_{0}^{1}d\alpha \,{\alpha }^{\gamma -2}\,{s}_{{\ell }}^{2}\left[\alpha ,e,{e}^{{\prime} }\right].\end{array}\end{eqnarray} \tag{ 124 }$$

In order to shorten the notation, let us introduce the dimensionless integrals

$$\begin{eqnarray}\begin{array}{rcl}{I}_{{\ell }}^{(1)}[p,e] & = & \displaystyle \int {{de}}^{{\prime} }\,{f}_{e}({e}^{{\prime} }){\displaystyle \int }_{0}^{1}d\alpha \,{\alpha }^{p}\,{s}_{{\ell }}^{2}\left[\alpha ,e,{e}^{{\prime} }\right],\\ {J}_{{\ell }}^{(1)}[p,e] & = & \displaystyle \int {{de}}^{{\prime} }\,{f}_{e}({e}^{{\prime} }){\displaystyle \int }_{0}^{1}d\alpha \,{\alpha }^{p}\,{s}_{{\ell }}^{2}\left[\alpha ,{e}^{{\prime} },e\right],\end{array}\end{eqnarray} \tag{ 125 }$$

where one should pay attention to the order of the arguments of ${s}_{{\ell }}^{2}$. This allows us then to rewrite Equation (122) as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}^{2}({\boldsymbol{K}})={{\rm{\Gamma }}}_{0}^{2}\,\displaystyle \frac{{f}_{{{\rm{\Gamma }}}^{2}}(e)}{1-{e}^{2}}{\left(\displaystyle \frac{a}{{a}_{0}}\right)}^{-\gamma },\end{eqnarray} \tag{ 126 }$$

where we introduce the amplitude Γ₀ and the dimensionless function ${f}_{{{\rm{\Gamma }}}^{2}}(e)$ as

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}_{0}^{2} & = & \displaystyle \frac{1}{4\pi }\,\displaystyle \frac{{{GN}}_{0}\langle {m}^{2}\rangle }{{M}_{\bullet }{a}_{0}^{3}},\\ {f}_{{{\rm{\Gamma }}}^{2}}(e) & = & \displaystyle \sum _{{\ell }}{B}_{{\ell }}\left\{{I}_{{\ell }}^{(1)}[\gamma -2,e]+{J}_{{\ell }}^{(1)}[2-\gamma ,e]\right\}.\end{array}\end{eqnarray} \tag{ 127 }$$

In Figure 10, we illustrate the dependence of ${f}_{{{\rm{\Gamma }}}^{2}}(e)$, assuming a thermal eccentricity distribution, f_e(e) = 2e, and different cusp profiles. We also note that Equation (127) is essentially identical to the result already presented in Equation (D5) of Kocsis & Tremaine (2015).

Zoom In Zoom Out Reset image size

Let us now pursue a similar approach to compute the coherence time, ${T}_{{\rm{c}}}^{{\rm{t}}}$, as introduced in Equation (39). When expanding the rhs of that equation, one gets

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Gamma }}}^{2}({\boldsymbol{K}})\,{T}_{{\rm{c}}}^{{\rm{t}}}({\boldsymbol{K}}) & = & {{\rm{\Gamma }}}_{0}\,\displaystyle \frac{{a}_{0}^{2}}{a(1-{e}^{2})}\\ & & \times \ \displaystyle \sum _{{\ell }}\displaystyle \frac{{B}_{{\ell }}}{\sqrt{{A}_{{\ell }}/2}}\displaystyle \int {{de}}^{{\prime} }\,\displaystyle \frac{\sqrt{1-{e}^{{\prime} 2}}\,{f}_{e}({e}^{{\prime} })}{\sqrt{{f}_{{{\rm{\Gamma }}}^{2}}({e}^{{\prime} })}}\\ & & \times \ \displaystyle \int {{da}}^{{\prime} }\,\displaystyle \frac{1}{{a}_{\mathrm{out}}^{2}}{\left(\displaystyle \frac{{a}^{{\prime} }}{{a}_{0}}\right)}^{\tfrac{4-\gamma }{2}}{s}_{{\ell }}^{2}\left[\tfrac{{a}_{\mathrm{in}}}{{a}_{\mathrm{out}}},{e}_{\mathrm{in}},{e}_{\mathrm{out}}\right].\end{array}\end{eqnarray} \tag{ 128 }$$

The integral over da' can be carried out following the same approach as in Equations (123) and (124) to get

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{0}^{a}{{da}}^{{\prime} }\,\displaystyle \frac{1}{{a}_{\mathrm{out}}^{2}}{\left(\displaystyle \frac{{a}^{{\prime} }}{{a}_{0}}\right)}^{\tfrac{4-\gamma }{2}}{s}_{{\ell }}^{2}\left[\tfrac{{a}_{\mathrm{in}}}{{a}_{\mathrm{out}}},{e}_{\mathrm{in}},{e}_{\mathrm{out}}\right]\\ \quad \ =\ \displaystyle \frac{1}{{a}_{0}}{\left(\displaystyle \frac{a}{{a}_{0}}\right)}^{\tfrac{2-\gamma }{2}}{\displaystyle \int }_{0}^{1}d\alpha \,{\alpha }^{\tfrac{4-\gamma }{2}}\,{s}_{{\ell }}^{2}\left[\alpha ,{e}^{{\prime} },e\right],\\ {\displaystyle \int }_{a}^{+\infty }{{da}}^{{\prime} }\,\displaystyle \frac{1}{{a}_{\mathrm{out}}^{2}}{\left(\displaystyle \frac{{a}^{{\prime} }}{{a}_{0}}\right)}^{\tfrac{4-\gamma }{2}}{s}_{{\ell }}^{2}\left[\tfrac{{a}_{\mathrm{in}}}{{a}_{\mathrm{out}}},{e}_{\mathrm{in}},{e}_{\mathrm{out}}\right]\\ \quad \ =\ \displaystyle \frac{1}{{a}_{0}}{\left(\displaystyle \frac{a}{{a}_{0}}\right)}^{\tfrac{2-\gamma }{2}}{\displaystyle \int }_{0}^{1}d\alpha \,{\alpha }^{\tfrac{\gamma -4}{2}}\,{s}_{{\ell }}^{2}\left[\alpha ,e,{e}^{{\prime} }\right].\end{array}\end{eqnarray} \tag{ 129 }$$

Similarly to Equation (125), in order to shorten the notation, we define the dimensionless integrals

$$\begin{eqnarray}\begin{array}{rcl}{I}_{{\ell }}^{(2)}[p,e] & = & \displaystyle \int {{de}}^{{\prime} }\,\displaystyle \frac{\sqrt{1-{e}^{{\prime} 2}}\,{f}_{e}({e}^{{\prime} })}{\sqrt{{f}_{{{\rm{\Gamma }}}^{2}}({e}^{{\prime} })}}{\displaystyle \int }_{0}^{1}d\alpha \,{\alpha }^{p}\,{s}_{{\ell }}^{2}\left[\alpha ,e,{e}^{{\prime} }\right],\\ {J}_{{\ell }}^{(2)}[p,e] & = & \displaystyle \int {{de}}^{{\prime} }\,\displaystyle \frac{\sqrt{1-{e}^{{\prime} 2}}\,{f}_{e}({e}^{{\prime} })}{\sqrt{{f}_{{{\rm{\Gamma }}}^{2}}({e}^{{\prime} })}}{\displaystyle \int }_{0}^{1}d\alpha \,{\alpha }^{p}\,{s}_{{\ell }}^{2}\left[\alpha ,{e}^{{\prime} },e\right],\end{array}\end{eqnarray} \tag{ 130 }$$

where once again, one should pay attention to the order of the arguments of ${s}_{{\ell }}^{2}$. Gathering all these elements, Equation (128) gives us the required expression for ${T}_{{\rm{c}}}^{{\rm{t}}}$. It reads

$$\begin{eqnarray}&&{T}_{{\rm{c}}}^{{\rm{t}}}({\boldsymbol{K}})={T}_{0}\,{f}_{T}(e){\left(\displaystyle \frac{a}{{a}_{0}}\right)}^{\gamma /2},\end{eqnarray} \tag{ 131 }$$

where we introduce the amplitude T₀ and the dimensionless function f_T (e) as

$$\begin{eqnarray}\begin{array}{rcl}{T}_{0} & = & \displaystyle \frac{1}{{{\rm{\Gamma }}}_{0}};\,{f}_{T}(e)=\displaystyle \frac{1}{{f}_{{{\rm{\Gamma }}}^{2}}(e)}\\ & & \times \ \displaystyle \sum _{{\ell }}\displaystyle \frac{{B}_{{\ell }}}{\sqrt{{A}_{{\ell }}/2}}\left\{{I}_{{\ell }}^{(2)}\left[\tfrac{\gamma -4}{2},e\right]+{J}_{{\ell }}^{(2)}\left[\tfrac{4-\gamma }{2},e\right]\right\}.\end{array}\end{eqnarray} \tag{ 132 }$$

In Figure 10, for a thermal eccentricity distribution, we illustrate how one can assume ${f}_{{{\rm{\Gamma }}}^{2}}^{1/2}(e)\,{f}_{T}(e)\simeq 0.4$ independent of e and the considered cups's power index γ.