Vector Resonant Relaxation of Stars around a Massive Black Hole

TL;DR

This paper develops an analytical framework for vector resonant relaxation of stars near a massive black hole, comparing it with simulations and providing a stochastic model for orbital reorientation.

Contribution

It introduces a first-principles derivation of potential fluctuation correlations and a stochastic scheme for orbital orientation evolution, advancing understanding of stellar dynamics near black holes.

Findings

Analytical correlation functions for potential fluctuations

Agreement between theory and N-body simulations

Quantitative estimates for Milky Way-like clusters

Abstract

In the vicinity of a massive black hole, stars move on precessing Keplerian orbits. The mutual stochastic gravitational torques between the stellar orbits drive a rapid reorientation of their orbital planes, through a process called vector resonant relaxation. We derive, from first principles, the correlation of the potential fluctuations in such a system, and the statistical properties of random walks undergone by the stellar orbital orientations. We compare this new analytical approach with effective $N$-body simulations. We also provide a simple scheme to generate the random walk of a test star's orbital orientation using a stochastic equation of motion. We finally present quantitative estimations of this process for a nuclear stellar cluster such as the one of the Milky Way.

Equations

$$H$$

$$=-L(\mathbf{K})\sum _{\begin{array}{c}\ell \ge 2\\ \mathrm{even}\end{array}}\sum _{m=-\ell }^{\ell }{M}_{\ell m}(\mathbf{K},t)Y_{\ell m}(\hat{\mathbf{L}}).$$

$$M_{\ell m}(\mathbf{K},t)=\sum_{i=1}^{N}\mathcal{J}_{\ell}\big{[}\mathbf{K},\mathbf{K}_{i}\big{]}\,Y_{\ell m}(\widehat{\mathbf{L}}_{i}(t)),$$

$$\frac{\partial \varphi }{\partial t}=\frac{\partial H}{\partial L_z};\frac{\partial L_z}{\partial t}=-\frac{\partial H}{\partial \varphi },$$

$$\frac{\mathrm{d}\mathbf{L}}{\mathrm{d}t}=\frac{\partial }{\partial \mathbf{L}}\times (H\mathbf{L})=L\sum _{\ell ,m}{M}_{\ell m}(\mathbf{K},t){\mathbf{X}}_{\ell m}(\hat{\mathbf{L}}(t)),$$

$$\phi (\hat{\mathbf{L}},\mathbf{K},t)=\sum _i{\delta }_{\mathrm{D}}(\hat{\mathbf{L}}-{\hat{\mathbf{L}}}_i(t)){\delta }_{\mathrm{D}}(\mathbf{K}-{\mathbf{K}}_i),$$

$$\frac{\partial \phi (\hat{\mathbf{L}},\mathbf{K},t)}{\partial t}=-$$

$$\times Y_{\ell m}({\hat{\mathbf{L}}}^{\prime }){\mathbf{X}}_{\ell m}(\hat{\mathbf{L}})\cdot \frac{\partial \phi (\hat{\mathbf{L}},\mathbf{K},t)}{\partial \hat{\mathbf{L}}},$$

$$\phi (\hat{\mathbf{L}},\mathbf{K},t)={\phi }_{\alpha }(\mathbf{K},t)Y_{\alpha }(\hat{\mathbf{L}}),$$

$${\phi }_{\alpha }(\mathbf{K},t)=\sum _i{\delta }_{\mathrm{D}}(\mathbf{K}-{\mathbf{K}}_i)Y_{\alpha }({\hat{\mathbf{L}}}_i(t)).$$

$$\frac{\partial {\phi }_{\alpha }(\mathbf{K},t)}{\partial t}=-\int \mathrm{d}{\mathbf{K}}^{\prime }{Q}_{\alpha \gamma \delta }(\mathbf{K},{\mathbf{K}}^{\prime }){\phi }_{\gamma }({\mathbf{K}}^{\prime },t){\phi }_{\delta }(\mathbf{K},t),$$

$$E_{\alpha \gamma \delta }=\int \mathrm{d}\hat{\mathbf{L}}{Y}_{\alpha }(\hat{\mathbf{L}}){\mathbf{X}}_{\gamma }(\hat{\mathbf{L}})\cdot \frac{\partial Y_{\delta }(\hat{\mathbf{L}})}{\partial \hat{\mathbf{L}}}.$$

$$C_{\alpha \beta }(\mathbf{K},{\mathbf{K}}^{\prime },t-t^{\prime })\equiv ⟨{\phi }_{\alpha }(\mathbf{K},t){\phi }_{\beta }({\mathbf{K}}^{\prime },t^{\prime })⟩,$$

$$C_{\alpha \beta }(\mathbf{K},{\mathbf{K}}^{\prime },0)={\delta }_{\alpha }^{\beta }{\delta }_{\mathrm{D}}(\mathbf{K}-{\mathbf{K}}^{\prime })n(\mathbf{K})$$

$$\displaystyle\frac{\partial^{2}}{\partial t^{2}}C_{\alpha\beta}(\mathbf{K},\mathbf{K}^{\prime},t)\bigg{|}_{t=0}\!\!\!=-\delta_{\alpha}^{\beta}\,A_{\ell_{\alpha}}\delta_{\mathrm{D}}(\mathbf{K}-\mathbf{K}^{\prime})\,n(\mathbf{K})\,\Gamma^{2}(\mathbf{K}),$$

$$A_{\ell }=\ell (\ell +1).$$

$${\Gamma }^2(\mathbf{K})=\int \mathrm{d}{\mathbf{K}}^{\prime }n({\mathbf{K}}^{\prime })\sum _{\ell }{B}_{\ell }{\mathcal{J}}_{\ell }^2[\mathbf{K},{\mathbf{K}}^{\prime }],$$

$$B_{\ell }=\frac{\ell (\ell +1)(2\ell +1)}{8\pi }.$$

$$C_{\alpha \beta }(\mathbf{K},{\mathbf{K}}^{\prime },t)={\delta }_{\alpha }^{\beta }{\delta }_{\mathrm{D}}(\mathbf{K}-{\mathbf{K}}^{\prime })C_{{\ell }_{\alpha }}(\mathbf{K},t),$$

$$C_{\ell }(\mathbf{K},t)=n(\mathbf{K}){\mathrm{e}}^{-\frac{A_{\ell }}{2}{(t/T_{\mathrm{c}}(\mathbf{K}))}^2},$$

$$T_{\mathrm{c}}(\mathbf{K})=\frac{1}{\Gamma (\mathbf{K})}.$$

$$C_{\alpha \beta }^{\mathrm{r}}(\mathbf{K},{\mathbf{K}}^{\prime },t-t^{\prime })\equiv {⟨{\phi }_{\alpha }(\mathbf{K},t){\phi }_{\beta }({\mathbf{K}}^{\prime },t^{\prime })⟩}_T,$$

$$⟨f{⟩}_T\equiv \frac{1}{T}{\int }_0^T\mathrm{d}tf(t)$$

$$C_{\alpha \beta }^{\mathrm{r}}(\mathbf{K},{\mathbf{K}}^{\prime },t)={\delta }_{\alpha }^{\beta }{\delta }_{\mathrm{D}}(\mathbf{K}-{\mathbf{K}}^{\prime })C_{{\ell }_{\alpha }}^{\mathrm{r}}(\mathbf{K},t),$$

$$C_{\ell }^{\mathrm{r}}(\mathbf{K},t)=n_{\ell }^{\mathrm{r}}(\mathbf{K}){\mathrm{e}}^{-\frac{A_{\ell }}{2}{(t/T_{\mathrm{c}}(\mathbf{K}))}^2}.$$

$$n_{{\ell }_{\alpha }}^{\mathrm{r}}(\mathbf{K})=\frac{1}{2{\ell }_{\alpha }+1}\sum _{m_{\alpha },\beta }\int \mathrm{d}{\mathbf{K}}_{\beta }{C}_{\alpha \beta }^{\mathrm{r}}({\mathbf{K}}_{\alpha },{\mathbf{K}}_{\beta },0),$$

$$\displaystyle\big{\langle}n_{\ell}^{\mathrm{r}}(\mathbf{K})\big{\rangle}$$

$${\phi }^{\mathrm{t}}(\hat{\mathbf{L}},t)={\delta }_{\mathrm{D}}(\hat{\mathbf{L}}-{\hat{\mathbf{L}}}_{\mathrm{t}}(t)),$$

$${\phi }_{\alpha }^{\mathrm{t}}(t)=Y_{\alpha }({\hat{\mathbf{L}}}_{\mathrm{t}}(t)).$$

$$\frac{\partial {\phi }_{\alpha }^{\mathrm{t}}(t)}{\partial t}$$

$$=-Q_{\alpha \delta }^{\mathrm{t}}(t){\phi }_{\delta }^{\mathrm{t}}(t),$$

$$\varphi^{\mathrm{t}}_{\alpha}(t)={\big{[}\mathrm{e}^{\Omega(t^{\prime},t)}\big{]}}_{\alpha\delta}\,\varphi^{\mathrm{t}}_{\delta}(t^{\prime}),$$

Full text

Galaxy: center - Galaxy: nucleus - galaxies: nuclei - gravitation - celestial mechanics

1 Introduction

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 Fig. 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 orientation of the stars get reshuffled, without changing the magnitude of their angular momentum nor their Keplerian energy, through a process called vector resonant relaxation (VRR) (Kocsis & Tremaine, 2015, and references therein), which is the focus on this work. Resonant torques’ coupling between the precessing stars then lead 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 the stars’ 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 binaries mergers (Hamers et al., 2018). Since its first presentation by Rauch & Tremaine (1996), VRR was studied numerically both with full (Eilon et al., 2009) and effective orbit-averaged (Kocsis & Tremaine, 2015) $N$-body simulations. More recently, Roupas et al. (2017); Takács & Kocsis (2018); Szölgyén & Kocsis (2018) studied the thermodynamical equilibria of VRR, and Fouvry et al. (2018) its axisymmetric limit.

In the present 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 then 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 effective 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. Finally, in Section 6, we use this new formalism to present the timescales associated with VRR in a nuclear stellar cluster similar to the Milky Way’s. We conclude in Section 7.

2 Model

We consider a set of $N$ stars orbiting a 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 semi-major 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_{\mathrm{p}}\!=\!(1-e)a}$ and ${r_{\mathrm{a}}\!=\!(1+e)a}$ with surface density ${\Sigma(r)\!=\!{[2\pi^{2}a\sqrt{r\!-\!r_{\mathrm{p}}}\sqrt{r_{\mathrm{a}}\!-\!r}]}^{-1}}$, where the semi-major axis $a$ and eccentricity $e$ can be assumed to be constant in time (see an illustration in Fig. 1).

Following this double orbit-average, one can associate to each star a set of conserved quantities ${\mathbf{K}=(m,a,e)}$ and a time-dependent normal vector ${\widehat{\mathbf{L}}}$, with $\bm{L}=L\widehat{\mathbf{L}}$ the orbital angular momentum. We introduce the spherical coordinates as ${(\theta,\phi)}$, so that ${\widehat{\mathbf{L}}=(\sqrt{1-u^{2}}\cos\phi,\sqrt{1-u^{2}}\sin\phi,u)}$, with ${u=\cos(\theta)}$. Studying VRR amounts then to studying the long-term dynamics of each star’s normal vector $\widehat{\mathbf{L}}$.

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

[TABLE]

with ${\mathbf{r}(t)}$, ${\mathbf{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 $\mathbf{K}$ and $\mathbf{K}_{i}$ their respective conserved parameters. In the second line of Eq. (1), we introduced the magnetizations

[TABLE]

where the coupling coefficients ${\mathcal{J}_{\ell}[\mathbf{K},\mathbf{K}_{j}]}$ are defined in Eq. (A1), and we used real spherical harmonics ${Y_{\ell m}(\widehat{\mathbf{L}})}$ (defined in Eq. (B1)). It is also important to note that only even harmonics with ${\ell\geq 2}$ contribute the particles’ dynamics, in virtue of the symmetries of the interaction.

Hamilton’s equations of motion read

[TABLE]

where ${L_{z}=Lu}$ is an action and $\phi$ its conjugated angle. The evolution of the angular momentum of a single test particle is given by

[TABLE]

where ${\mathbf{X}_{\ell m}(\widehat{\mathbf{L}})=\widehat{\mathbf{L}}\!\times\!\partial Y_{\ell m}(\widehat{\mathbf{L}})/\partial\widehat{\mathbf{L}}}$ are the real vector spherical harmonics. Because $L$ is constant, one has ${\mathrm{d}\widehat{\mathbf{L}}/\mathrm{d}t\!=\!L^{-1}\mathrm{d}\bm{L}/\mathrm{d}t}$.

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

[TABLE]

with${\delta_{\mathrm{D}}}$ the Dirac delta, ${\delta_{\mathrm{D}}(\widehat{\mathbf{L}}\!-\!\widehat{\mathbf{L}}_{i})\!=\!\delta_{\mathrm{D}}(u\!-\!u_{i})\delta_{\mathrm{D}}(\phi\!-\!\phi_{i})}$, and ${\delta_{\mathrm{D}}(\mathbf{K}\!-\!\mathbf{K}_{i})\!=\!\delta_{\mathrm{D}}(m\!-\!m_{i})\delta_{\mathrm{D}}(a\!-\!a_{i})\delta_{\mathrm{D}}(e\!-\!e_{i})}$, with the associated volumes ${\mathrm{d}\widehat{\mathbf{L}}=\mathrm{d}u\mathrm{d}\phi}$ and ${\mathrm{d}\mathbf{K}=\mathrm{d}m\mathrm{d}a\mathrm{d}e}$. The continuity equation, ${\partial\varphi/\partial t=-\partial/\partial\widehat{\mathbf{L}}\!\cdot\![\varphi\,\partial\widehat{\mathbf{L}}/\partial t]}$, gives us then

[TABLE]

where we used the fact that the vector spherical harmonics satisfy ${\partial/\partial\widehat{\mathbf{L}}\!\cdot\!\mathbf{X}(\widehat{\mathbf{L}})=0}$. This equation can subsequently be developed in spherical harmonics, by writing

[TABLE]

where the sum over the index ${\alpha=(\ell_{\alpha},m_{\alpha})}$ is implied, and we introduced

[TABLE]

When expanded in spherical harmonics, Eq. (6) becomes

[TABLE]

where the sums over the harmonic indices ${(\gamma,\delta)}$ are implied, and we introduced the time-independent coupling tensor ${Q_{\alpha\gamma\delta}(\mathbf{K},\mathbf{K}^{\prime})\!=\!\mathcal{J}_{\ell_{\gamma}}[\mathbf{K},\mathbf{K}^{\prime}]\,E_{\alpha\gamma\delta}}$, with $E_{\alpha\gamma\delta}$ the (real) Elsasser coefficients (James, 1973) (see Appendix B for their properties)

[TABLE]

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}(\mathbf{K},t)}$.

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

3 The correlation function of the noise

As a first step towards the characterization of the correlated stochastic dynamics of one star in that system, we focus our interest 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 Eq. (22)) for the time dependences of this same correlation function.

The harmonic coefficients ${\varphi_{\alpha}(\mathbf{K},t)}$ (Eq. 8) describe the full state of the ${N\gg 1}$ particles system 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

[TABLE]

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 and generically decreases to zero on a timescale larger than some coherence time $T_{\mathrm{c}}$. As a result, as a first approximation, it is therefore reasonable to replace ${C_{\alpha\beta}(\mathbf{K},\mathbf{K}^{\prime},t-t^{\prime})}$ by a Gaussian function, tailored to match the function’s behavior for ${t\ll T_{\mathrm{c}}}$, see Fig. 3 for a justification.

In Appendix D, we compute the first two derivatives of the correlation function, and we show in Eqs. (D4) and (D8) that

[TABLE]

and

[TABLE]

with the coefficient

[TABLE]

In Eq. (13), we introduced ${n(\mathbf{K})}$, the DF of the stars’ ${\mathbf{K}\!=\!(a,e,m)}$ parameters, that satisfies the normalization convention ${\!\int\!\mathrm{d}\widehat{\mathbf{L}}\mathrm{d}\mathbf{K}n(\mathbf{K})\!=\!N}$. We also introduced the decay rate of the correlation function $\Gamma(\mathbf{K})$ as

[TABLE]

with the coefficient

[TABLE]

Gathering Eqs. (12) and (13), we can approximate the correlation function $C_{\alpha\beta}(\mathbf{K},\mathbf{K}^{\prime},t)$ by

[TABLE]

where ${C_{\ell}(\mathbf{K},t)}$ is a function decaying like a Gaussian

[TABLE]

where we introduced the torque time

[TABLE]

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 Fig. 3 we compare this estimate to the correlations measured in the $N$-body 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 $N$-body 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 Fig. 3. Here, such a calculation is made intricate by our accounting of the $\ell$- and $\mathbf{K}$-dependence of the pairwise coupling, as embodied by the sum over $\ell$ and the integral over ${\mathrm{d}\mathbf{K}^{\prime}}$ in Eq. (9).

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 ‘$\mathrm{r}$’, we therefore define the time-averaged correlation function

[TABLE]

where

[TABLE]

stands for the time average over some long timescale $T$. As previously, we will assume that $C_{\alpha\beta}^{\mathrm{r}}(\mathbf{K},\mathbf{K}^{\prime},t)$ can be replaced by

[TABLE]

with the Gaussian time dependence

[TABLE]

Here, we defined an effective isotropic amplitude ${n_{\ell}^{\mathrm{r}}(\mathbf{K})}$ as the mean value over $m_{\alpha}$, so that

[TABLE]

and from Eqs. (12), we have

[TABLE]

that is independent of the considered harmonic. It is important to note that ${n_{\ell}^{\mathrm{r}}(\mathbf{K})}$ varies between different realizations. In Eq. (F23), we illustrate how one can compute its variance, and show how this originates from the constraint of total energy conservation.

In Eq. (23), we assumed, for simplicity, that the torque time, ${T_{\mathrm{c}}(\mathbf{K})}$, is independent of the considered realization. These various choices ensure that the ansatz from Eq. (22) satisfies the constraints from Eqs. (12) and (13) when ensemble-averaged. As highlighted by Eqs. (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_{\mathrm{D}}(\mathbf{K}-\mathbf{K}^{\prime})}$). Following Eq. (23), we note that the time-dependence of this correlation is controlled by both an isotropic amplitude, ${n_{\ell}^{\mathrm{r}}(\mathbf{K})}$, and a torque time, ${T_{\mathrm{c}}(\mathbf{K})}$, that depend on the considered parameter $\mathbf{K}$.

4 The random walk of a test particle

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 Eq. (17), our goal is now to investigate the stochastic dynamics of one given test particle embedded in that fluctuating environment. In Fig. 4, we illustrate one such random walk by highlighting the time evolution of the orientation of a single particle in one fiducial simulation.

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 $\mathbf{K}_{\mathrm{t}}$, and its orientation at time $t$ with ${\widehat{\mathbf{L}}_{\mathrm{t}}(t)}$. Similarly to Eq. (5), the current orientation of the test particle is fully characterized by the single particle DF

[TABLE]

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

[TABLE]

Since ${\varphi^{\mathrm{t}}_{(1,m)}(\widehat{\mathbf{L}}_{\mathrm{t}})\propto(y_{\mathrm{t}},z_{\mathrm{t}},x_{\mathrm{t}})\propto\widehat{\mathbf{L}}_{\mathrm{t}}}$, characterizing the random walk of the test particle on the sphere as in Fig. 4, requires the knowledge of the correlation properties of ${\varphi^{\mathrm{t}}_{\alpha}(t)}$ for ${\ell_{\alpha}=1}$.

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

[TABLE]

where ${\varphi_{\gamma}(\mathbf{K},t)}$ is the harmonic coefficients of the background particles’ field, as defined in Eq. (8), and we introduced ${Q^{\mathrm{t}}_{\alpha\delta}(t)\!=\!\!\int\!\mathrm{d}\mathbf{K}\,Q_{\alpha\gamma\delta}(\mathbf{K}_{\mathrm{t}},\mathbf{K})\,\varphi_{\gamma}(\mathbf{K},t)}$, that 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}(\mathbf{K},t)}$), whose correlation properties were investigated in the previous section. We also emphasize that this forcing term also depends on $\mathbf{K}_{\mathrm{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^{\mathrm{t}}_{\alpha}(t)$. In order to guarantee well-behaved asymptotics of the test particle’s motion for large times, we approach the resolution of Eq. (28) via 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

[TABLE]

where the matrix ${\Omega(t^{\prime},t)}$ is constructed as a series expansion of the form ${\Omega(t^{\prime},t)=\sum_{k\geq 1}\Omega_{k}(t^{\prime},t)}$, whose first terms are

[TABLE]

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

[TABLE]

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

[TABLE]

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

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

[TABLE]

where we introduced the dimensionless function

[TABLE]

As can be seen from the time-dependence of the exponent in equation (33), one can note that on short timescales, ${t\ll T_{\mathrm{c}}}$, the correlation $C^{\mathrm{t},\mathrm{r}}_{\alpha\beta}$ decays like a Gaussian and the motion of the particle is ballistic (e.g., ${(\Delta\widehat{\mathbf{L}})^{2}\propto t^{2}}$). On that short timescales, the random walk of the test particle is analogous to the one induced by a time-independent fluctuation. On long timescales, ${t\gg T_{\mathrm{c}}}$, the correlation decays exponentially in time and the motion of the test star is diffusive (e.g., ${(\Delta\widehat{\mathbf{L}})^{2}\propto t}$). On these long timescales, the random walk of the test particle is analogous to the one induced by fluctuations $\delta_{\mathrm{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 $\mathbf{K}$, Eq. (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 the ones of Eq. (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

[TABLE]

where the Gaussian noise ${\bm{\eta}(t)}$ is a ${3D}$ vector of zero mean, ${\langle\eta_{i}(t)\rangle=0}$, and follows ${\langle\eta_{i}(t)\,\eta_{j}(t^{\prime})\rangle=\delta_{ij}\,\mathrm{e}^{-{[(t-t^{\prime})/T_{\mathrm{c}}^{\mathrm{t}}]}^{2}}}$. We will then choose the amplitude $\Gamma_{\mathrm{t}}$ and coherence time $T_{\mathrm{c}}^{\mathrm{t}}$ by matching the short and long timescales behavior of the test particle’s correlation function with the ballistic and diffusive regimes of the generic result from Eq. (33), an approach already used in Hamers et al. (2018).

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

[TABLE]

with ${A_{\ell}=\ell(\ell+1)}$. By matching the ballistic and diffusive regimes of Eq. (4) with the ones of Eq. (33), we may then constrain the amplitude, $\Gamma_{\mathrm{t}}$, and coherence time, $T_{\mathrm{c}}^{\mathrm{t}}$, of the toy model of Eq. (35).

The amplitude, ${\Gamma_{\mathrm{t}}(\mathbf{K}_{\mathrm{t}})}$, varies from realization to realization and is given by

[TABLE]

When ensemble-averaged over realizations, this amplitude becomes

[TABLE]

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

[TABLE]

where for simplicity we assumed, similarly to Eq. (23), that the coherence time, ${T_{\mathrm{c}}^{\mathrm{t}}(\mathbf{K}_{\mathrm{t}})}$, is independent of the considered realization.

Equation (4) 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 the fluctuating noise from the background particles. The test particle’s random walk is characterized by the two quantities ${(\Gamma_{\mathrm{t}}^{2},T_{\mathrm{c}}^{\mathrm{t}})}$, that both depend on the test particle’s parameters $\mathbf{K}_{\mathrm{t}}$. The coefficient ${\Gamma_{\mathrm{t}}}$ controls the amplitude of the test particle’s initial ballistic motion, while the coherence time ${T_{\mathrm{c}}^{\mathrm{t}}}$ controls the timescale after which the test particle enters the diffusive regime. One strength of the present formalism is that, following Eqs. (37) and (39), one now has at one’s disposal explicit expressions for these two parameters. These coefficients can then easily be computed for various cluster models (by varying the DF ${n(\mathbf{K})}$) and various test particles (by varying $\mathbf{K}_{\mathrm{t}}$). Considering the same test particle as in Fig. 4, we illustrate in Fig. 5 one random walk generated using the Langevin equation (35).

However, as we had already emphasized in Eq. (22), it is important to note that the present system suffers from being non-ergodic, i.e. ensemble averages and time averages cannot be interverted. This is highlighted by the fact that even for the exact same test particle (i.e. same $\mathbf{K}_{\mathrm{t}}$), the amplitude ${\Gamma_{\mathrm{t}}}$ varies from realization to realization. In Appendix F, we compute the associated variance, as given by Eq. (F25), 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 Eq. (4). This is illustrated in Fig. 6, and to shorten the main text, we detailed in Appendix H.2 the procedure followed to obtain that figure.

In that figure, we note that the $N$-body measurements and the analytical prediction from Eq. (4) agree both on short timescales but also on timescales longer than the coherence time $T_{W_{\mathrm{t}}}^{\mathrm{t}}$ (defined in Eq. (H11)). As already stressed in Eq. (33), the second panel of Fig. 6 clearly exhibits the two successive regimes of evolution, namely ballistic for ${t\ll T_{\mathrm{c}}^{\mathrm{t}}}$ and diffusive for ${t\gg T_{\mathrm{c}}^{\mathrm{t}}}$. This same panel also emphasizes the importance of accounting for the variance in ${\Gamma_{\mathrm{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^{\mathrm{t}}_{(1,m)}(t)\propto\widehat{\mathbf{L}}_{\mathrm{t}}}$, Fig. 6 also offers then an illustration of the behavior of ${t\mapsto\langle\widehat{\mathbf{L}}_{\mathrm{t}}(t)\,\widehat{\mathbf{L}}_{\mathrm{t}}(0)\rangle}$. It is also straightforward to adapt that prediction to different test stars (by changing $\mathbf{K}_{\mathrm{t}}$) or different galactic nuclei (by changing the DF ${n(\mathbf{K})}$).

5 The self-consistency of the noise

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 Eq. (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 Eq. (20), and replace ${\varphi_{\alpha}(\mathbf{K},t)}$ by its definition in terms of a discrete sum over particles, as in Eq. (8), so that

[TABLE]

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 Eq. (40) becomes

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

and one finally gets the self-consistent differential equation111Similar self-consistent differential equations were first obtained in Taylor & McNamara (1971) (see Eq. (22) therein) in the context of plasma diffusion in the guiding center limit.

[TABLE]

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\!\mathrm{d}\mathbf{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 Eq. (18). This (motivated) ansatz may then be reinjected in the r.h.s. of the self-consistency relation from Eq. (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 Fig. 3, where we show how one can better match the late-time properties of the system’s noise through this iterative process.

6 Application

As an illustration of the present formalism, let us now consider the case of a stellar cusp distribution similar to the one of SgrA*. The mass of the MBH is taken to be ${M_{\bullet}\!=\!4\!\times\!10^{6}M_{\odot}}$, and for simplicity we consider a single-mass stellar population of individual mass ${m_{\star}\!=\!1M_{\odot}}$. 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_{0}\!=\!g(\gamma)N(\!<\!a_{0})}$, with

[TABLE]

and ${N(\!<\!a_{0})}$ the physical number of stars within a sphere of radius $a_{0}$ from the center. For the numerical application, we assume that ${a_{0}\!=\!r_{\mathrm{h}}\!=\!2\,\mathrm{pc}}$ and ${N(\!<\!a_{0})\!=\!4\!\times\!10^{6}}$. 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\pi)}$.

In Appendix I, we show that the amplitude ${\Gamma^{2}}$, defined in Eq. (15) and characterizing the ballistic regime of the orientation’s random walk, follows the power law distribution

[TABLE]

where ${P(a)=2\pi(a^{3}/(GM_{\bullet}))^{1/2}}$ is the orbital period, and ${f_{\Gamma^{2}}(e)\simeq 0.15}$ is a dimensionless eccentricity function defined in Eq. (I7) and illustrated in Fig. 10. In Fig. 7, we illustrate the dependence of the torque time, ${1/\Gamma}$, for circular orbits of different semi-major axes and for different cusp’s 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).

As shown in Fig. 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_{\mathrm{c}}^{\mathrm{t}}$, defined in Eq. (39) and characterizing the diffusive regime of the orientation’s random walk. It follows the power law distribution

[TABLE]

where the dimensionless eccentricity function ${f_{T}(e)}$ is given in Eq. (I12) and illustrated in Fig. 10. In that same figure, we note that for a thermal eccentricity distribution and a cusp’s power index ${1\leq\gamma\leq 2}$, one can assume that ${f_{\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

[TABLE]

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 Eq. (35), as one only has to estimate the test particle’s torque time ${1/\Gamma(a,e)}$, as the associated coherence time, ${T_{\mathrm{c}}^{\mathrm{t}}(a,e)}$, follows immediately.

7 Conclusion

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 of 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 interverted.

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 Eq. (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 Eq. (4), 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 timescales predictions. Throughout the text, all the predictions were compared with detailed effective $N$-body 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 the one of SgrA*.

The present paper is only a first step towards a complete theory of VRR, and we list below some possible tracks for future developments. 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 for quantitative predictions for the properties of these random walks, as obtained in Eq. (4). 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}[\mathbf{K},\mathbf{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. JBF 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. BB is supported by membership from Martin A. and Helen Chooljian at the Institute for Advanced Study.

Appendix A The coupling coefficients

Similarly to Eq. (10) in Kocsis & Tremaine (2015), we define the coupling coefficients222With this convention, the definition of $\mathcal{J}_{ij\ell}$ from Eq. (10) of Kocsis & Tremaine (2015) is recovered by ${\mathcal{J}_{ij\ell}=L[\mathbf{K}_{i}]\,\mathcal{J}_{\ell}[\mathbf{K}_{i},\mathbf{K}_{j}]\,(2\ell+1)/(4\pi)}$.

[TABLE]

where ${L[\mathbf{K}]=m\sqrt{GM_{\bullet}a(1-e^{2})}}$ is the magnitude of the angular momentum. We also introduced “$\mathrm{out}$” (resp. “$\mathrm{in}$”) as the index $i$ or $j$ with the larger (resp. smaller) semi-major axis, and defined accordingly the ratio ${\alpha=a_{\mathrm{in}}/a_{\mathrm{out}}\leq 1}$. With these notations, the dimensionless coefficients ${s_{\ell}[\alpha,e_{\mathrm{in}},e_{\mathrm{out}}]}$ are given by

[TABLE]

with ${P_{\ell}(u)}$ the usual Legendre polynomials. Because they are independent of the details of the considered system, the coefficients ${s_{\ell}[\alpha,e_{\mathrm{in}},e_{\mathrm{out}}]}$ can be precomputed on a grid to hasten the numerical evaluation of ${\mathcal{J}_{\ell}[\mathbf{K}_{i},\mathbf{K}_{j}]}$. For our fiducial simulations, these coefficients were pre-computed on a linear ${3D}$ grid in ${(\alpha,e_{\mathrm{in}},e_{\mathrm{out}})}$ consisting of $200^{3}$ elements, with ${10^{-2}\leq\alpha\leq 1}$ and ${0\leq e_{\mathrm{in}},e_{\mathrm{out}}\leq 0.99}$. We refer to Fig. 1 in Kocsis & Tremaine (2015) for an illustration of the behavior of these coefficients.

Appendix B The Elsasser coefficients

In this Appendix, we follow James (1973); Ivers & Phillips (2008) and detail some of the properties of the Elsasser 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

[TABLE]

with ${P_{\ell}^{m}(u)}$ the usual associated Legendre polynomials, and the coefficients {K_{\ell}^{m}=\big{[}\tfrac{2\ell+1}{4\pi}\frac{(\ell-m)!}{(\ell+m)!}\big{]}^{1/2}}. With this convention, the spherical harmonics follow the normalization ${\!\int\!\mathrm{d}\widehat{\mathbf{L}}\,Y_{\ell m}Y_{\ell^{\prime}m^{\prime}}=\delta_{\ell}^{\ell^{\prime}}\delta_{m}^{m^{\prime}}}$.

Following Eq. (5.5) of Ivers & Phillips (2008), the real Elsasser coefficients, as defined in Eq. (10), can be decomposed as

[TABLE]

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

[TABLE]

where we introduced the Wigner ${3j}$-symbols (Arfken et al., 2005), and defined

[TABLE]

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

[TABLE]

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

[TABLE]

The Elsasser coefficients satisfy various exclusion rules (James, 1973). In particular, for $E_{\alpha\gamma\delta}$ to be non-zero, one has to satisfy

[TABLE]

These coefficients also follow the symmetry relations ${E_{\alpha\delta\gamma}=E_{\gamma\alpha\delta}=-E_{\alpha\gamma\delta}}$.

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

[TABLE]

Appendix C $N$-body simulations

In this Appendix, we briefly detail the effective $N$-body simulations to which our analytical results are compared. To simulate a system of $N$ interacting particles, the starting point is the evolution Eq. (4), that can be used for each of the $N$ particles. In that form, we note that the velocity vector, ${\mathrm{d}\widehat{\mathbf{L}}_{i}/\mathrm{d}t}$, is expressed only as a function of the current location of the particle, ${\widehat{\mathbf{L}}_{i}(t)}$, and the instantaneous particle’s magnetizations, ${M_{\ell m}(\mathbf{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_{\rm max}^{2})}$, with $\ell_{\rm max}$ the maximum harmonic number considered in the pairwise interaction333 We note the similarities between the evolution Eq. (4) and the evolution equation for the Hamiltonian Mean Field (HMF) model (Antoni & Ruffo, 1995). The main differences are that (i) the present Hamiltonian has no kinetic term, (ii) the magnetizations depend on two harmonic indices ${(\ell,m)}$ (versus one for the HMF model), (iii) the magnetizations are not global but vary from one particle to another because of the parameters $\mathbf{K}$. As a result, the complexity of the integration of the HMF model scales like ${\mathcal{O}(N\ell_{\rm max})}$, versus ${\mathcal{O}(N^{2}\ell_{\rm max}^{2})}$ for the present model. . 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.

Our $N$-body 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 Eq. (2), (iii) computing the velocity fields in Eq. (4), (iv) advancing all the particles’ orientation for one timestep. The real spherical harmonics are computed using a reccurence relation for the renormalized associated Legendre polynomials (see Eq. (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\phi)}$) 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 ${\mathrm{d}\widehat{\mathbf{L}}_{i}/\mathrm{d}t}$ are computed, particles are advanced for a timestep $h$, using a fourth-order Runge-Kutta integrator (see Eq. (17.1.3) Press et al., 2007).

All the derivations presented in the main text are illustrated with comparisons with this direct $N$-body approach. We consider a system composed of ${N=10^{3}}$ stars, and assume that the particles’ conserved quantities ${\mathbf{K}_{i}=(m_{i},a_{i},e_{i})}$ satisfy ${m=m_{\mathrm{min}}}$, ${a_{\mathrm{min}}\leq a\leq a_{\mathrm{max}}}$, and ${e_{\mathrm{min}}\leq e\leq e_{\mathrm{max}}}$. Our units are chosen so that ${m_{\mathrm{min}}=a_{\mathrm{min}}=G=1}$, and we pick ${a_{\mathrm{max}}/a_{\mathrm{min}}=100}$, ${e_{\mathrm{min}}=0}$, and ${e_{\mathrm{max}}=0.3}$. These parameters are drawn independently one from another, according to probability distribution functions proportional to ${(\delta_{\mathrm{D}}(m-m_{\mathrm{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 ${\ell_{\rm 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}=|\mathrm{d}\widehat{\mathbf{L}}_{i}/\mathrm{d}t|=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\!\text{Min}_{i}[t_{\tau}^{i}]}$. With such choices, integrating the system for one timestep takes approximately ${1\mathrm{s}}$ on a single core, and simulations are carried out for ${2\!\times\!10^{5}}$ timesteps.

Appendix D Computing the derivatives of the noise correlation

D.1 Computing ensemble averages

We are generically interested in computing ensemble averages at the initial time of the form ${\langle\varphi_{\alpha}(\mathbf{K},0)\,\varphi_{\beta}(\mathbf{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{\mathbf{L}})=1/(4\pi)}$, while we assume that their parameter $\mathbf{K}$ is drawn according to a PDF, ${g(\mathbf{K})}$, normalized so that ${\!\int\!\mathrm{d}\mathbf{K}\,g(\mathbf{K})=1}$.

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

[TABLE]

where non-zero terms only come from ${i=j}$, and we introduced the connected average as

[TABLE]

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

[TABLE]

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 Eq. (11).

For the value at the initial time, we can write

[TABLE]

where we followed Eq. (D1) to compute the last average, and introduced ${n(\mathbf{K})=g(\mathbf{K})\,N/(4\pi)}$ as the DF of the stars’ parameters satisfying the normalization ${\!\int\!\mathrm{d}\widehat{\mathbf{L}}\mathrm{d}\mathbf{K}\,n(\mathbf{K})\!=\!N}$. We emphasize that to compute Eq. (D4), 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\!\sim\!\langle\varphi_{\gamma}\,\varphi_{\delta}\varphi_{\beta}\rangle\!=\!0}$, where we used the quadratic evolution Eq. (9) once. As it involves an odd number of fields, this correlation is equal to zero, as imposed by Eq. (D3).

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

[TABLE]

where we injected the evolution Eq. (9) twice. As shown in Eq. (D3), in the limit of Gaussian fluctuations, the average term can be computed by keeping only averages of pairs. As ${E_{\alpha\gamma\gamma}=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

[TABLE]

where we used ${E_{\beta\delta\gamma}=-E_{\beta\gamma\delta}}$, and introduced

[TABLE]

Following Eq. (B8), one can now perform the sums over $m_{\gamma}$ and $m_{\delta}$ in Eq. (D6). The VRR interactions being limited to even harmonic numbers $\ell$, we may impose at this stage that $\ell_{\alpha}$ is even. Glancing back at the constraint {C2} from Eq. (B7), we note that ${\ell_{\alpha}\!+\!\ell_{\gamma}\!+\!\ell_{\delta}}$ has to be odd, so that the term ${\Lambda_{\gamma\delta}(\mathbf{K},\mathbf{K}^{\prime})}$ never contributes to Eq. (D6). For $\ell_{\alpha}$ even, Eq. (D6) becomes

[TABLE]

where the sum over $\ell_{\delta}$ was performed following Eq. (B8), and the decay rate ${\Gamma^{2}(\mathbf{K})}$ is given by Eq. (15).

Appendix E Computing the properties of the random walk

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

[TABLE]

To compute the time average of ${\mathrm{e}^{\Omega(t)}}$, we rely on the cumulant theorem,

[TABLE]

where ${\mu_{n}=\left\langle\Omega^{n}\right\rangle_{T}}$ are the moment matrices, and $\kappa_{n}$ are the cumulant matrices, with the first two ones given by ${\kappa_{1}=\mu_{1}}$ and ${\kappa_{2}=\mu_{2}-\mu_{1}^{2}}$. Here, we compute the time average of ${\mathrm{e}^{\Omega(t)}}$ by keeping only terms that are at most second order in $Q^{\mathrm{t}}$, so that only ${\langle\Omega_{1}(t)\rangle_{T}\!\propto\!Q^{\mathrm{t}}}$ and ${\langle\Omega_{1}^{2}(t)\rangle_{T}\!\propto\!\langle\Omega_{2}(t)\rangle_{T}\!\propto\!{(Q^{\mathrm{t}})}^{2}}$ contribute. We note that since ${\langle\varphi_{\alpha}(\mathbf{K},t)\rangle_{T}\!=\!0}$ and, from stationarity, ${\langle Q^{\mathrm{t}}(t)\,Q^{\mathrm{t}}(t^{\prime})\rangle_{T}\!=\!\langle Q^{\mathrm{t}}(t^{\prime})\,Q^{\mathrm{t}}(t)\rangle_{T}}$, both $\langle\Omega_{1}(t)\rangle_{T}$ and $\langle\Omega_{2}(t)\rangle_{T}$ vanish. As a result, at the order considered here, only the second cumulant $\kappa_{2}=\langle\Omega_{1}^{2}(t)\rangle_{T}$ is non-zero, and from the cumulant theorem we obtain

[TABLE]

Gathering Eqs. (E1) and (E3) together, we can write the test particle’s correlation function as

[TABLE]

Using Eq. (30), we can write

[TABLE]

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 Eq. (22). Using the sum identities from Eq. (B8), one gets

[TABLE]

where ${\chi_{\ell}^{\mathrm{r}}(\mathbf{K},t)}$ is a double time integral of the noise correlation

[TABLE]

with the dimensionless function ${\chi(\tau)}$ defined in Eq. (4). Since the matrix $\langle\Omega_{1}^{2}(t)\rangle_{T}$ is diagonal, one can straightforwardly compute its exponential, as required by Eq. (E3). This allows us to finally recast the correlation of the test particle’s random motion from Eq. (E4) under the form of Eq. (33).

Appendix F Computing the variance of the noise amplitude

In this Appendix, we compute the ensemble-averaged variance of the amplitude of the density fluctuations, ${C_{\alpha\beta}^{\mathrm{r}}(\mathbf{K}_{\alpha},\mathbf{K}_{\beta},0)}$, as introduced in Eq. (20). Our goal is to compute an expression of the form

[TABLE]

Because only even harmonics contribute to the interactions, we can limit ourselves to ${2\leq\ell_{\alpha},\ell_{\gamma}}$ even. As estimated in Eq. (23), we know that for ${|t-t^{\prime}|\gg T_{\mathrm{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 interverting time- and ensemble-averages. Provided that these constraints are satisfied, in the two-dimensional time integral from Eq. (F1), one can note that the particles are uncorrelated between $t$ and $t^{\prime}$ on a surface of size ${(T-T_{\mathrm{c}})^{2}}$, while they are correlated on a surface of size ${TT_{\mathrm{c}}}$. As a result, as long as ${T\gg T_{\mathrm{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 Eq. (F1), in the presence of these constraints. At time $t$, the state of the system is fully characterized by the set of all fields ${\bm{\varphi}\!=\!\{\varphi_{\alpha}(\mathbf{K}_{\alpha},t)\}}$, and similarly, at time $t^{\prime}$, the state of the system is fully characterized by ${\bm{\varphi}^{\prime}\!=\!\{\varphi^{\prime}_{\alpha}(\mathbf{K}_{\alpha},t)\}}$. We may then use $\bm{\varphi}$ and $\bm{\varphi}^{\prime}$ as the random variables over which averages are carried out. Following Eq. (12), we have

[TABLE]

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

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

[TABLE]

with ${F_{\bm{\theta}}(\bm{\theta})}$ the PDF of the constraints $\bm{\theta}$. Therefore, we can write

[TABLE]

In that view, Eq. (F1) can be recast as

[TABLE]

where we introduced

[TABLE]

with ${\langle\,\cdot\,\rangle}$ standing for the ensemble average where the fields $\varphi$ are drawn according to the Gaussian statistics of ${F(\bm{\varphi})}$. Conveniently, in that form, Eq. (F5) allows us to carry out independently the averages over $\bm{\varphi}$ and $\bm{\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 ${(\bm{\theta}_{0})}$, the conservation of the system’s total angular momentum ${(\bm{\theta}_{1})}$, and the conservation of the system’s total energy ${(\bm{\theta}_{2})}$. Luckily, these can all be expressed as simple functions of the fields $\bm{\varphi}$. They read

[TABLE]

with ${L[\mathbf{K}]\!=\!m\sqrt{GM_{\bullet}a(1-e^{2})}}$ the norm of the angular momentum and ${H_{\ell}[\mathbf{K},\mathbf{K}^{\prime}]=L[\mathbf{K}]\,\mathcal{J}_{\ell}[\mathbf{K},\mathbf{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 ${\ell=0}$ for the conservation of $\mathbf{K}$, ${\ell=1}$ for the angular momentum, and ${2\leq\ell}$ 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 Eq. (F5), as we will now argue.

Since only even harmonics contribute to the interactions, we can restrict ourselves to ${2\leq\ell_{\alpha}}$ even when computing ${\mathbb{F}_{\alpha\beta}(\mathbf{K}_{\alpha},\mathbf{K}_{\beta},\bm{\theta})}$. If ${\ell_{\beta}=0,1}$, Eq. (F6) can be rewritten as

[TABLE]

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(\bm{\varphi}_{\ell\geq 2})}$, is an even function of the fields, the last bracket in Eq. (F8) is equal to zero. As a result, we can assume that ${\ell_{\beta}\geq 2}$. In that case, Eq. (F6) becomes

[TABLE]

where we introduced

[TABLE]

As a result, for ${\ell_{\alpha},\ell_{\beta},\ell_{\gamma},\ell_{\delta}\geq 2}$, this allows us to rewrite the needed correlation from Eq. (F5) as

[TABLE]

where we got rid of all occurences of the constraints $\bm{\theta}_{0}$ and $\bm{\theta}_{1}$, using the fact that their PDFs satisfy ${\!\int\!\mathrm{d}\bm{\theta}_{0}\,F_{\bm{\theta}_{0}}\!=\!\!\int\!\mathrm{d}\bm{\theta}_{1}F_{\bm{\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 Eq. (F10). In the present Gaussian limit, we can rely on Novikov theorem (Novikov, 1965) to compute it.444We do not repeat here the general theory of Novikov theorem (Novikov, 1965). We refer to Hänggi (1978) for non-Gaussian noises, to Garcia-Ojalvo & Sancho (1999) for spatially extended noises, and to Fouvry & Bar-Or (2018) for an example of application in stellar dynamics. One gets

[TABLE]

where the first cumulant is absent because ${\ell_{\alpha}\geq 2}$, so that ${\langle\varphi_{\alpha}(\mathbf{K}_{\alpha})\rangle=0}$, and only the second cumulant remains as the fields are assumed to be Gaussian. In Eq. (F12), the sum (resp. integral) over $\mu$ (resp. ${\mathrm{d}\mathbf{K}_{\mu}}$) runs over all the fields. The functional gradient appearing in the last term can be computed as

[TABLE]

where we used the fundamental relation ${\delta\varphi_{\beta}(\mathbf{K}_{\beta})/\delta\varphi_{\mu}(\mathbf{K}_{\mu})=\delta_{\mu}^{\beta}\delta_{\mathrm{D}}(\mathbf{K}_{\beta}-\mathbf{K}_{\mu})}$. Glancing back at the definition of the energy in Eq. (F7), we can also write

[TABLE]

Injecting these results into Eq. (F12) and using the Gaussian statistics from Eq. (F2), we obtain a self-consistent integro-differential equation for ${\mathbb{F}_{\alpha\beta}(\mathbf{K}_{\alpha},\mathbf{K}_{\beta},E)}$, namely

[TABLE]

where we used that ${\langle\delta_{\mathrm{D}}(E(\bm{\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_{\ell}\to\varepsilon H_{\ell}}$ in Eq. (F15), and consider the expansion

[TABLE]

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

[TABLE]

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

[TABLE]

where, for simplicity, we did not repeat the arguments ${(\mathbf{K}_{\alpha},\mathbf{K}_{\beta})}$ and ${(\mathbf{K}_{\gamma},\mathbf{K}_{\delta})}$. The zeroth-order term is straightforward to compute, and gives

[TABLE]

using ${\!\int\!\mathrm{d}E\,F_{E}(E)=1}$. It is straighforward to show that the first-order term satisfies

[TABLE]

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

[TABLE]

where ${\Delta E}$ is obtained after straightforward manipulations of the energy integrals, and reads

[TABLE]

It is important to note that ${\Delta E}$ is a single number that depends only on the total energy PDF, ${F_{E}(E)}$, and therefore only on the considered DF  ${n(\mathbf{K})}$. In Appendix G, we detail how the needed integral from Eq. (F22) can be estimated. Equation (F21) 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 Eq. (24), we can get the variance of ${n_{\ell}^{r}(\mathbf{K})}$. It reads

[TABLE]

where we introduced the dimensionless function ${M_{\ell}(\mathbf{K})}$ as

[TABLE]

The final step of this Appendix is to compute the variance of $\Gamma_{\mathrm{t}}^{2}$, as defined in Eq. (37). We get

[TABLE]

with ${\Delta\Gamma_{\mathrm{t}}^{2}\!=\!\Gamma_{\mathrm{t}}^{2}\!-\!\Gamma^{2}}$. Equation (F25) is the final result of this section. It expresses the variance (over realizations) of ${\Gamma_{\mathrm{t}}^{2}}$, the amplitude of the random walk of a given test star. It is important to note that this since ${(\Gamma^{2})^{2}}$ and ${\langle(\Delta\Gamma_{\mathrm{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 Fig. 6, let us illustrate the effect associated with this non-zero variance of ${\Gamma_{\mathrm{t}}^{2}}$ in our fiducial simulations. We consider the same window, ${W_{\mathrm{t}}(\mathbf{K}_{\mathrm{t}})}$, as in Eq. (H8). For each test particle falling in that window, we measure the correlation function ${\langle\varphi_{\alpha}^{\mathrm{t}}(t)\,\varphi_{\alpha}^{\mathrm{t}}(0)\rangle_{T}}$ (e.g. for ${\ell_{\alpha}=1}$). The second-order time derivative at ${t=0}$ of this correlation function is directly proportional to ${\Gamma_{\mathrm{t},\mathrm{num.}}^{2}}$ (see Eq. (4)), that we can therefore measure numerically. In Fig. 8, we represent the distribution of these numerically measured initial values, ${\Gamma_{\mathrm{t},\mathrm{num.}}^{2}}$, and illustrate how these amplitudes vary from realizations to realizations.

To capture this variance effect seen in Fig. 8, we may use our estimation of the variance of $\Gamma_{\mathrm{t}}^{2}$ obtained in Eq. (F25). To do so, for every test particle falling in the window, we compute ${(\Gamma^{2},\langle(\Delta\Gamma_{\mathrm{t}}^{2})^{2}\rangle)}$, following Eqs. (38) and (F25). Having determined this mean and variance, one can draw a sample of ${\Gamma_{\mathrm{t}}^{2}}$ according to a PDF that shares the same first two cumulants. Similarly to Eq. (G5), 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 Fig. 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 ${\Gamma_{\mathrm{t}}^{2}}$.

Appendix G Computing the variance of the energy distribution

In this Appendix, we briefly detail how one can estimate the energy spread ${\Delta E}$ introduced in Eq. (F22). Our convention for the definition of the total energy is spelled out in Eq. (F7). Using the two-point statistics from Eq. (D2), the expected mean value of the energy reads

[TABLE]

where we defined

[TABLE]

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

[TABLE]

To compute the average term appearing in the r.h.s., we follow Eq. (D3), placing ourselves in the limit of Gaussian random fields so that only connected averages involving two fields remain. We get

[TABLE]

where we used the symmetry relation ${H_{\ell}[\mathbf{K},\mathbf{K}^{\prime}]=H_{\ell}[\mathbf{K}^{\prime},\mathbf{K}]}$, and introduced

[TABLE]

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

[TABLE]

with ${\kappa=\mu^{2}/\sigma^{2}}$. For our fiducial numerical system, we find ${\langle E\rangle\simeq 6.1\!\times\!10^{-3}}$ and ${\kappa\simeq 81}$. In Fig. 9, we illustrate the statistical distribution of the system’s energy, as well as the approximation from Eq. (G6).

Finally, for a chi-squared PDF as in Eq. (G6), one can explicitly compute ${\Delta E}$, as defined in Eq. (F22), to get

[TABLE]

We note that this integral is well-behaved only for ${\kappa>2}$. This is an artefact coming from the perturbative expansion introduced in Eq. (F16). For our fiducial model, we find ${\Delta E\simeq 6.7\!\times\!10^{-4}}$.

Appendix H Computing averages over window

In this Appendix, we briefly detail the procedures used in Figs. 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 Fig. 3 used to illustrate Eq. (22). One of the key insight from this equation is that to any particle (of parameter $\mathbf{K}$), we can associate the pair ${(n(\mathbf{K}),T_{\mathrm{c}}(\mathbf{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 $\mathbf{K}$-averaged fields ${\overline{\varphi}_{\alpha}(t)\!=\!\!\int\!\mathrm{d}\mathbf{K}\,W(\mathbf{K})\,\varphi_{\alpha}(\mathbf{K},t)}$, with ${W(\mathbf{K})}$ a window function defined555As detailed in Appendix C, our fiducial simulations are single-mass, so that ${n(\mathbf{K})\!\propto\!\delta_{\mathrm{D}}(m-m_{\mathrm{min}})\,g_{a}(a)\,g_{e}(e)}$. As a consequence, for the definition of the window function ${W(\mathbf{K})}$ in Eq. (H1) to be meaningful, we do not account for the Dirac delta in mass present in ${n(\mathbf{K})}$. as

[TABLE]

with ${(C_{\mathrm{\min}},T_{\mathrm{c}}^{\mathrm{\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{\mathbf{L}}_{i}(t))}$, with the sum limited to the particles with ${(n(\mathbf{K}),T_{\mathrm{c}}(\mathbf{K}))}$ in the vicinity of ${(C_{\mathrm{min}},T_{\mathrm{c}}^{\mathrm{min}})}$. Introducing ${C_{\ell_{\alpha},W}^{\mathrm{r}}(t-t^{\prime})\equiv\langle\overline{\varphi}_{\alpha}(t)\,\overline{\varphi}_{\alpha}(t^{\prime})\rangle_{T}}$, Eq. (22) immediately gives

[TABLE]

where the function ${C_{\ell}^{\mathrm{r}}(\mathbf{K},t)}$ follows the Gaussian ansatz from Eq. (23). When averaged over realizations, Eq. (H2) can be approximated with the Gaussian dependence

[TABLE]

where we introduced the amplitude ${C_{W}}$ and torque time ${T_{W}}$ as

[TABLE]

Equation (H3) is the analytical Gaussian prediction illustrated in Fig. 3.

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

[TABLE]

with $C_{W}$ as in Eq. (H4), and where the amplitude $\Gamma_{W}^{2}$ and coherence time ${T_{W}^{\mathrm{t}}}$ are given by

[TABLE]

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

[TABLE]

H.2 Correlation of the random walks

Let us briefly detail the method followed to obtain Fig. 6, used to illustrate the result from Eq. (4). Following the independence hypothesis from Eq. (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 one particular particle.

One important insight from Eq. (4) is that to any test particle (of parameter $\mathbf{K}_{\mathrm{t}}$), we can associate the pair ${(\Gamma^{2}(\mathbf{K}_{\mathrm{t}}),T_{\mathrm{c}}^{\mathrm{t}}(\mathbf{K}_{t}))}$, that characterizes the correlation properties of its random walk in orientation. Similarly to Eq. (H1), in order to investigate these random walks, it is convenient to introduce the window function ${W_{\mathrm{t}}(\mathbf{K}_{\mathrm{t}})}$ as

[TABLE]

where ${(\Gamma_{\mathrm{min}}^{2},T_{\mathrm{t}}^{\mathrm{min}})}$ are the typical amplitude and coherence time of the considered test particles, and $\varepsilon_{W_{\mathrm{t}}}$ is a dimensionless parameter controlling the size of the window. We may then define the window-averaged correlation function (that can easily be measured in the numerical simulations) as

[TABLE]

where ${\langle\,\cdot\,\rangle_{W_{\mathrm{t}}}}$ stands for the mean over the test particles of a given realization falling in the window ${W_{\mathrm{t}}(\mathbf{K}_{\mathrm{t}})}$, similarly to Eq. (H7). Following Eq. (4), we also added a prefactor ${4\pi}$ to ensure that this correlation is between [math] and $1$.

If one does not account for the variance in ${\Gamma_{\mathrm{t}}^{2}}$ (see Eq. (F25)), a first (naive) prediction for Eq. (H9) can be obtained from Eq. (4) by restricting ourselves only to the ensemble-averaged mean prediction-s. This gives

[TABLE]

where the amplitude, ${\Gamma_{W_{\mathrm{t}}}^{2}}$, and coherence time, ${T_{W_{\mathrm{t}}}}$, are computed by direct averages over the particles falling in the window ${W_{\mathrm{t}}(\mathbf{K}_{\mathrm{t}})}$, so that

[TABLE]

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

[TABLE]

Here, the amplitude and coherence time from Eq. (H11) become

[TABLE]

where ${\mathrm{PDF}[\mu,\sigma^{2}]}$ returns a sample from a chosen PDF of mean $\mu$ and variance $\sigma^{2}$ (which we chose to be a chi-squared PDF as in Fig. 8). Both predictions from Eqs. (H10) and (H12) are illustrated in Fig. 6, for the particular harmonic ${\ell=1}$.

Appendix I The case of a power law distribution

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

[TABLE]

where we recall that “$\mathrm{out}$” (resp. “$\mathrm{in}$”) labels the star with the larger (resp. smaller) semi-major axis, and we introduced the dimensionless ratio ${\alpha=a_{\mathrm{in}}/a_{\mathrm{out}}}$.

Following Eq. (15), we can then compute the amplitude, $\Gamma^{2}$, of the background fluctuations. It reads

[TABLE]

where we introduced the second moment of the mass distribution ${\langle m^{2}\rangle\!=\!\!\int\!\mathrm{d}mf_{m}(m)\,m^{2}}$. The integral over $a^{\prime}$ in Eq. (I2) can then be split into two regions, ${a^{\prime}\leq a}$ and ${a^{\prime}\geq a}$. For the first region, we write

[TABLE]

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

[TABLE]

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

[TABLE]

where one should pay attention to the order of the arguments of ${s_{\ell}^{2}}$. This allows us then to rewrite Eq. (I2) as

[TABLE]

where we introduced the amplitude $\Gamma_{0}$ and the dimensionless function ${f_{\Gamma^{2}}(e)}$ as

[TABLE]

In Fig. 10, we illustrate the dependence of ${f_{\Gamma^{2}}(e)}$, assuming a thermal eccentricity distribution, ${f_{e}(e)=2e}$, and different cusp’s profiles.

Let us now pursue a similar approach to compute the coherence time, $T_{\mathrm{c}}^{\mathrm{t}}$, as introduced in Eq. (39). When expanding the r.h.s. of that equation, one gets

[TABLE]

The integral over ${\mathrm{d}a^{\prime}}$ can be carried out following the same approach as in Eqs. (I3) and (I4) to get

[TABLE]

Similarly to Eq. (I5), in order to shorten the notations, we define the dimensionless integrals

[TABLE]

where once again, one should pay attention to the order of the arguments of $s_{\ell}^{2}$. Gathering all these elements, Eq. (I8) gives us the needed expression of $T_{\mathrm{c}}^{\mathrm{t}}$. It reads

[TABLE]

where we introduced the amplitude $T_{0}$ and the dimensionless function ${f_{T}(e)}$ as

[TABLE]

In Fig. 10, for a thermal eccentricity distribution, we illustrate how one can assume ${f_{\Gamma^{2}}^{1/2}(e)\,f_{T}(e)\simeq 0.4}$ independent of $e$ and the considered cups’s power index $\gamma$.