© The Authors 2024

1. Introduction

Rotation is ubiquitous in stellar systems. In effect, it provides a source of free energy, allowing clusters to efficiently reshuffle their orbital structure towards more likely configurations. Yet, historically, the study of globular clusters has been mainly focused on isotropic, non-rotating, old globular clusters (Aarseth et al. 1974; Spitzer 1975; Cohn 1979; Trager et al. 1995; Miocchi et al. 2013). The reasons behind such simplifications are two-fold: (i) naturally, it is easier numerically and analytically to neglect the effect of rotation; (ii) spherical isotropic models – for example, the King models (King 1966) or the Wilson models (Wilson 1975) – provided a satisfactory zeroth-order description of the main observed dynamical properties of globular clusters (see, e.g., McLaughlin & van der Marel 2005).

The last decade has seen the extraction of new data – for instance, HST (Bellini et al. 2017) and Gaia DR2 (Bianchini 2018; Sollima et al. 2019). These surveys gave the astrophysical community access to numerous and detailed observations of the internal kinematics of several globular clusters of the Milky Way (Bianchini et al. 2013; Fabricius et al. 2014; Watkins et al. 2015; Ferraro et al. 2018; Kamann et al. 2018), as well as a quantification of the degree of velocity anisotropy (Jindal et al. 2019). Using these new data sets, the historical highly symmetric cluster models are not satisfactory anymore.

Therefore, a secular theory that describes the evolution of globular clusters and accounts for their rotation is needed to describe their long-term evolution. In particular, the Fokker-Planck theory was used to probe the impact of rotation of the long-term evolution of such clusters, both before (Einsel & Spurzem 1999) and after (Kim et al. 2002) core collapse. Furthermore, observations show that the angular momentum distribution measured in Galactic clusters retain the signature of their formation process (Lanzoni et al. 2018). While N-body simulations are able to reproduce these results (Tiongco et al. 2016, 2022), the historical context and the complexity of the problem have led to few analytical explorations (Geyer et al. 1983; White & Shawl 1987; Kontizas et al. 1989). These remain scarce even now (see, e.g., Stetson et al. 2019; Rozier et al. 2019; Livernois et al. 2022).

In isolated systems, a non-zero total angular momentum, that is, the presence of rotation, can have a significant impact on the cluster’s long-term evolution, for instance, through the “gravo-gyro catastrophe” (Hachisu 1979; Ernst et al. 2007). This phenomenon has been observed in a range of rotating systems, including gas cylinders (Inagaki & Hachisu 1978; Hachisu 1979), gaseous discs (Hachisu 1982), flattened (quasi-spherical) star clusters (Akiyama & Sugimoto 1989; Einsel & Spurzem 1999; Ernst et al. 2007), and clusters with embedded black holes (BHs) (Fiestas & Spurzem 2010; Kamlah et al. 2022). Let us however note that many of these studies used a rotating King model to study the impact of rotation (see, e.g., Einsel & Spurzem 1999; Varri & Bertin 2012). Therein, changing the amount of rotation impacts the density profile: this makes comparisons amongst different models less clear.

The concurrent occurrence of internal rotation and a spectrum of stellar masses can result in the formation of an oblate core of fast rotating heavy masses (Kim et al. 2004; Tiongco et al. 2021). More precisely, the orbital inclinations of the heaviest stars align with respect to one another, inducing a mass segregation in the distribution of orbital inclinations (Szölgyén et al. 2019). The generally agreed explanation for this phenomenon is resonant relaxation and resonant friction (Rauch & Tremaine 1996; Meiron & Kocsis 2019). While this effect concerns globular clusters, nuclear cluster with a dominant massive BH (Szölgyén & Kocsis 2018; Foote et al. 2020; Gruzinov et al. 2020; Magnan et al. 2022; Ginat et al. 2023) also display this spontaneous alignment.

The objective of this paper is to explore the impact of rotation on the long-term evolution of spherically symmetric globular clusters. To achieve this goal, we extend the orbit average analysis of anisotropic Plummer globular clusters from Tep et al. (2022) to the case of rotating clusters. By performing tailored N-body simulations, we quantify the validity of Chandrasekhar’s non-resonant relaxation (NR) theory (Chandrasekhar 1943), and assess the importance of the gravo-gyro catastrophe as well as that of collective effects (i.e., the self-amplification of the star’s gravitational response, see Heyvaerts 2010).

This paper is organised as follows. In Section 2, we extend the NR theory to rotating clusters. In Section 3, we follow the long-term evolution of a series of rotating clusters using N-body simulations. We study the impact of rotation on core collapse and on the distribution of orbital inclinations. We compare the NR prediction with N-body simulations in Sections 4 (in-plane diffusion) and 5 (out-of-plane diffusion). Finally, we discuss our results in Section 6.

2. Non-resonant relaxation

We consider a self-gravitating globular cluster composed of N stars of individual mass m​ = ​M/N, with M the cluster’s total mass. We follow the cluster’s evolution through the total distribution function (DF) in (r, v) space, F_rot​ = ​F_rot(r, v), with r the position and v the velocity. All normalisations are taken to be the same as in Section 2 of Tep et al. (2022), hereafter T22. In addition to quasi-stationarity, we assume the cluster to be spherically symmetric – hence, with planar unperturbed orbits – and in rotation. Since we are interested in orbital distortion, it is convenient to monitor the cluster’s evolution in action space, via F_rot​ = ​F_rot(J), where J = (J_r, L, L_z) are the specific action coordinates. Here, J_r is the radial action, L the norm of the angular momentum vector and L_z its projection along the z-axis. We also introduce the orbital inclination through cos I​ = ​L_z/L. Because actions are integrals of motion, we use them to label the cluster’s orbits and track their deformation over time.

The long-term evolution of the cluster’s DF is described by the FP equation (see, e.g., section 7.4 of Binney & Tremaine 2008)

$$\begin{array}{cc}\hfill \frac{\partial F_{\mathrm{rot}}(\mathit{J},t)}{\partial t}& =-\frac{\partial }{\partial \mathit{J}}·\mathcal{F}(\mathit{J})\hfill \\ \hfill & =-\frac{\partial }{\partial \mathit{J}}·[{\mathit{D}}_1(\mathit{J})F_{\mathrm{rot}}(\mathit{J})-\frac{1}{2}\frac{\partial }{\partial \mathit{J}}·({\mathit{D}}_2(\mathit{J})F_{\mathrm{rot}}(\mathit{J}))],\hfill \end{array}$$(1)

with ℱ(J) the action space flux. Therein, we find the first-order diffusion coefficient, D₁, and the diffusion matrix, D₂. Both of them describe the distortion of orbits in action space. The next two sections briefly describe the main steps needed to compute this diffusion flux in Equation (1).

2.1. Local velocity deflection coefficients

Let us consider a test star of mass m and velocity v. As a result of the cluster’s finite number of constituents, this test star is subject to perturbations around its mean field trajectory, driving an irreversible diffusion of its velocity. We call this long-term relaxation process, sourced by successive, uncorrelated pairwise deflections, the NR relaxation. Following Appendix L of Binney & Tremaine (2008) (p. 836), the corresponding local velocity deflection coefficients generically read

$$\begin{array}{cc}\hfill ⟨\mathrm{\Delta }{v}_i⟩& =-2A{\displaystyle \int \mathrm{d}}{\mathit{v}}^{\prime }\frac{w_i}{w^3}{F}_{\mathrm{rot}}(\mathit{r},{\mathit{v}}^{\prime }),\hfill \end{array}$$(2a)

$$\begin{array}{cc}\hfill ⟨\mathrm{\Delta }{v}_i\mathrm{\Delta }{v}_j⟩& =A{\displaystyle \int \mathrm{d}}{\mathit{v}}^{\prime }\frac{w^2{\delta }_{\mathit{ij}}-w_i{w}_j}{w^3}{F}_{\mathrm{rot}}(\mathit{r},{\mathit{v}}^{\prime }),\hfill \end{array}$$(2b)

where 1​ ≤ ​i, j​ ≤ ​3 run over the three directions of the coordinate system. We also defined A​ = ​4πmG²​lnΛ and the relative velocity w​ = ​v − v′. Here, lnΛ stands for the Coulomb logarithm, set to Λ​ = ​0.11N as is usual for single-mass globular clusters (Giersz & Heggie 1994; Heggie & Hut 2003). In principle, the rewriting of these expressions using Rosenbluth potentials is feasible, as outlined by Rosenbluth et al. (1957). Such a program, in a non-rotating but anisotropic cluster, was pursued in T22.

In Appendix A, we improve upon this approach in the context of rotating clusters. In particular, we show how a rewriting rather based on Equations (2) improves the numerical stability and ensures the positivity of the second-order diffusion coefficients. More precisely, we write Equations (2) as

$$\begin{array}{c}\hfill \left[\begin{array}{c}⟨\mathrm{\Delta }{v}_{\Vert }⟩\\ ⟨{(\mathrm{\Delta }{v}_{\Vert })}^2⟩\\ ⟨{(\mathrm{\Delta }{v}_{\perp })}^2⟩\end{array}\right]=A{\displaystyle \int \mathrm{d}}w\mathrm{d}\vartheta \mathrm{d}\varphi sin\vartheta \left[\begin{array}{c}-2cos\vartheta \\ w{sin}^2\vartheta \\ w(1+{cos}^2\vartheta )\end{array}\right]F_{\mathrm{rot}},\end{array}$$(3)

where Δv_∥ (resp. Δv_⊥) is the local velocity deflection along (resp. perpendicular to) the test star’s trajectory and we have used polar coordinates (w, ϑ, ϕ) with the z-axis parallel to v (see Fig. A.1). Appendix A.1 also details the arguments at which the DF, $F_{\text{rot}}(\mathit{r},{\mathit{v}}^{\prime })\text{}=\text{}{F}_{\text{rot}}(E^{\prime },L^{\prime },L_z^{\prime })$, must be evaluated.

Equation (3) is one of the main results of the present work. In particular, Equation (3) possesses a few key advantages compared to the expressions given in Equations (3) of T22, which are recovered by the present approach: (i) no integrable singular denominators remain; (ii) no gradients of F_rot are required during the computation; (iii) the positivity of the second-order coefficients is ensured; (iv) this equation applies to a wider range of clusters. Ultimately, Equation (3) compactly accounts for all the two-body deflections from the cluster’s stars onto the test star.

2.2. Orbit average and secular evolution

Because we are interested in the diffusion of orbits, we must consider the effect of NR on orbital invariants, namely the energy E and the angular momenta L and L_z. In Appendix B, we expand upon T22 and detail how the local diffusion coefficients in (E, L, L_z) may be computed from the local velocity coefficients given in Equation (3).

Once the local diffusion coefficients are known, they can be averaged along the unperturbed mean field orbit of the test star. At this stage, our accounting of rotation adds some complexity. Indeed, the orbit average now involves an intricate two-dimensional integral spanning the radial range of the test orbit (r​ ∈ ​[r_p, r_a]) and its angular phase (φ​ ∈ ​[0, 2π]) within the orbital plane. This is illustrated in Fig. 1.

Fig. 1. Illustration of the orbit average in physical space. The test star (in black, part of whose orbit is shown in red) is averaged over all its available positions in its orbital plane, within the radial range of the orbit (r​ ∈ ​[rp, ra]) and its angular range (φ​ ∈ ​[0, 2π]), weighted by the surface density Ωr/(2π2|vr|). Here, the frame is that of the right panel of Fig. A.2.

In practice, computing the orbit average diffusion coefficients involves computing expressions of the form

$$\begin{array}{c}\hfill D_X=\frac{{\mathrm{\Omega }}_r}{\pi }{\displaystyle {\int }_{r_{\mathrm{p}}}^{r_{\mathrm{a}}}}\frac{\mathrm{d}r}{|v_{\mathrm{r}}|}{\displaystyle {\int }_0^{2\pi }}\frac{\mathrm{d}\phi }{2\pi }⟨\mathrm{\Delta }X⟩(r,\phi ),\end{array}$$(4)

where X runs over all the diffusion coefficients to compute, that is, runs over singlets and pairs from {E, L, L_z}. Here, Ω_r is the radial frequency and v_r the radial velocity that is evaluated along the mean field orbit. We refer to Appendix C for detailed expressions. In particular, in Appendix C.1, we show how the angular average can be performed explicitly for some classes of rotating models. And, in Appendix C.2, we use the same approach as in appendix F2 of T22 to compute the radial average using a stable numerical integration scheme.

The final stage of the calculation is to convert the orbit-averaged diffusion coefficients in (E, L, L_z) into ones in J​ = ​(J_r, L, L_z). This is rather straightforward, as detailed in Appendix D. Lastly, these diffusion coefficients are those used to evaluate the FP Equation (1). Along the same lines, one can similarly describe the dynamics within the coordinates J_c​ = ​(J_r, L, cos I).

3. Long-term relaxation

We wish to study the impact of rotation on the long-term evolution of rotating, spherically symmetric¹, anisotropic clusters. In this first section, we focus on direct N-body simulations to get some insight on the dynamics at play. We detail the corresponding numerical setup in Appendix E.

First, let us describe the classes of clusters considered. The in-plane distribution of orbits (i.e., when projecting in the (J_r, L) space) follows the same Plummer DFs (see, e.g., Dejonghe 1987) as in Section 3.1 of T22. Velocity anisotropy is encoded with the dimensionless parameter q, with q = 0 corresponding to an isotropic distribution, and q​ > ​0 (resp. q​ < ​0) associated with radially (resp. tangentially) anisotropic velocity distributions.

Let us introduce rotation in these models, while leaving the mean potential invariant. To do so, we follow the Lynden-Bell daemon (LBD) (Lynden-Bell 1960) and consider

$$\begin{array}{c}\hfill F_{\mathrm{rot}}(\mathit{J})=F_{\mathrm{tot}}(J_{\mathrm{r}},L)(1+\alpha \mathrm{sgn}[L_z/L]),\end{array}$$(5)

where F_tot(J_r, L) is the non-rotating DF and α a dimensionless parameter between 0 and 1. Physically, this parameter corresponds to converting a fraction α of retrograde orbits (L_z​ < ​0) into prograde orbits (L_z​ > ​0). Importantly, since L_z​ ↦ ​sgn[L_z/L] is an odd function, we stress that F_tot and F_rot generate the exact same potential.

When considering orbital inclinations via J_c​ = ​(J_r, L, cos I) (see Appendix D.3), the LBD yields the DF in J_c space²

$$\begin{array}{c}\hfill F({\mathit{J}}_{\mathrm{c}})=L{F}_{\mathrm{tot}}(J_{\mathrm{r}},L)(1+\alpha \mathrm{sgn}[cosI]).\end{array}$$(6)

Integrating this equation over cos I, one recovers the reduced DF in (J_r, L), F(J_r, L)​ = ​2LF_tot(J_r, L), used to describe non-rotating clusters (see, e.g., Hamilton et al. 2018).

3.1. Sphericity of the cluster

The NR theory presented in Section 2 assumes that the cluster remains spherically symmetric throughout. However, Rozier et al. (2019) showed that (sufficiently) rotating clusters can harbour unstable modes (see, e.g., Fig. 8 therein). To avoid this complication, we restricted ourselves to (linearly) stable rotating clusters, that is, clusters that would remain spherically symmetric.

To probe the conservation of spherical symmetry, we investigateed the “sphericity”, h, of clusters, as defined in Appendix F in terms of the eigenvalue ratio of the density-weighted inertia tensor. This is illustrated in Fig. 2.

Fig. 2. Sphericity, h, of a sample of rotating anisotropic clusters, as defined in Appendix F, from N-body simulations for N​ = ​105. Each case has been averaged over 50 realisations. We represent radially anisotropic clusters (q​ = ​1) in red, isotropic clusters (q​ = ​0) in yellow, and tangentially anisotropic ones (q​ = ​ − 6) in blue. For each anisotropy, we consider three rotating parameters α​ = ​0.1, 0.25, 0.5 (ordered from dark to light colors). In agreement with Rozier et al. (2019), some of these clusters are unstable – namely (q, α)​ = ​(1, 0.5), and in a smaller fashion (q, α)​ = ​(1, 0.25) – while the others are linearly stable and remain spherically symmetric.

Except for the radially-anisotropic and rapidly-rotating cluster (q, α)​ = ​(1, 0.5), all clusters remain approximately spherical. This is in agreement with the measurements from Rozier et al. (2019), as this particular cluster falls into the region of linear instability. Therefore, we shall focus our interest on parameters (q, α) for which the clusters remain spherically symmetric.

3.2. The gravo-gyro catastrophe

The typical relaxation timescale of a globular cluster can be estimated through the half-mass relaxation time (see, e.g., Section 14 of Heggie & Hut 2003, for more details), defined by

$$\begin{array}{c}\hfill t_{\mathrm{rh}}=\frac{0.138N{r}_{\mathrm{h}}^{3/2}}{{(GM)}^{1/2}ln(0.11N)},\end{array}$$(7)

with r_h the half-mass radius. For a cluster with N​ = ​10⁵ stars, this yields t_rh​ = ​994 HU, where HU stands for the Hénon units (Heggie & Mathieu 1986) in which G​ = ​M​ = ​R_v​ = ​1, with R_v the virial radius³. In practice, with N​ = ​10⁵ running a numerical simulation up to t​ ∼ ​t_rh took about one day of computation, see Appendix E. In Breen et al. (2017), core collapse occurs at t​ ∼ ​​17t_rh for an isotropic cluster (see Table 1 therein for the dependence with respect to velocity anisotropy). Hence, integrating clusters with N​ = ​10⁵ up to core collapse is not reasonably feasible. In this section we therefore scaled down the size of the clusters to N​ = ​10⁴ stars. In that case, t_rh​ = ​132 HU, and core collapse was numerically reached in about ten hours. Since we focused on measuring the core radius – a very integrated quantity – the quality of the numerical measurements was not too much degraded by this use of a smaller value for N. In practice, measurements were averaged over 50 runs.

Figure 3 represents the time evolution of the averaged core radius, R_c, defined as (see, e.g., Breen et al. 2017; T22)

$$\begin{array}{c}\hfill R_{\mathrm{c}}^2={\displaystyle \sum _i}{r}_i^2{\rho }_i^2/{\displaystyle \sum _i}{\rho }_i^2,\end{array}$$(8)

Fig. 3. Core radius as a function of time, as measured in N-body simulations, with trh defined in equation (7). Each case has been averaged over 50 realisations with N​ = ​104. Increasing the rotation strength α slightly reduces the time of core collapse. Nevertheless, the impact of rotation (i.e. the gravo-gyro catastrophe) is surely not as pronounced as what was observed in, e.g., Einsel & Spurzem (1999).

where r_i is the radial position of star i and ρ_i an estimator of the density at r_i, as defined in NBODY6++GPU (Wang et al. 2015). For the usual sets of anisotropies and rotations (q​ = ​1, 0, −6 and α​ = ​0, 0.1, 0.25, 0.5).

Interestingly, for these rotating clusters, the impact of rotation is definitely not as important as the one reported in Einsel & Spurzem (1999) in rotating King models (see Figure 2 therein). Appendix G reproduces the main parameters characterising these models. For example, in Einsel & Spurzem (1999), assuming W₀​ = ​6.0, a non-rotating King cluster collapses at t ∼ 12t_rh, while a rotating one, with ω₀/Ω₀ = 0.4, can collapse as early as t ∼ 9t_rh. Here, for Plummer spheres, we do not find any such stark impact of rotation.

Nonetheless, Figure 3 does not contradict the measurements from Einsel & Spurzem (1999). Indeed, the parameters ω₀ (King) and α (Plummer) do not parametrise rotation in the same way. As Einsel & Spurzem (1999) varies ω₀, the mean density profile gets modified. Indeed, the rotating DF of the King model (Equation G.1) cannot be decomposed into a fixed, rotation-free even part, and an odd part in L_z (see, e.g., Dejonghe 1986). As such, we argue that the LBD approach ensures a fairer comparison between models to isolate the distinctive impact of rotation.

Of course, to better assess the universality of these observations, it would be worthwhile to perform similar experiments with other parametrisations for the rotation (see Appendix K). This is left for future investigations.

3.3. In-plane vs out-of-plane diffusion

Using once again N-body simulations with N​ = ​10⁴ stars, let us finally investigate the typical features of the in-plane diffusion (i.e. relaxation in (J_r, L)) and out-of-plane diffusion (i.e. relaxation in cos I).

First, Fig. 4 reports on the time evolution of the DF in L, that is, the in-plane relaxation of the distribution of angular momenta. There, we recover the imprints of core collapse visible through the slow overall contraction of the distribution. But importantly, we clearly note that rotation only (very) weakly impacts this relaxation.

Fig. 4. Evolution of the angular momentum DF, F(L), as measured in N-body simulations with N​ = ​104, with trh defined in Equation (7). For each parameter, the latest time is close to the time of core collapse. Each panel is averaged over 50 realisations. Here, F(L) redistributes towards lower angular momenta, with details depending on the initial velocity anisotropy. Rotation does not impact strongly the in-plane relaxation.

We can now compare this relaxation with the out-of-plane one. Figure 5 reports on the time evolution of the DF in cos I. As could have been expected, the clusters tend towards smoother distribution of inclinations, meaning that relaxation tends to erase discontinuities. The sharp discontinuity at cos I​ = ​0 has already been washed away after 1 t_rh, regardless of the initial amplitude of the discontinuity. We refer to Appendix K.2 for a comparison of the early relaxation between a discontinuous and a smooth distribution of orbital inclinations. Note that the redistribution of prograde and retrograde orbits observed in Fig. 5 does not conflict with the conservation of the cluster’s total angular momentum. Indeed, the average of cos I​ = ​L_z/L is not constrained by any global invariance, contrary to L_z. As such, even if the number of particles with positive and negative cos I (hence L_z) changes, the conservation of the total angular momentum is ensured by a modulation of the norm of each star’s angular momentum vector.

Fig. 5. Same as Fig. 4 but for the DF in orbital inclinations, F(cos I). The out-of-plane relaxation does not depend much on the initial velocity anisotropy. As time evolves, the discontinuity of the LBD (Equation (5)) at cos I​ = ​0 is (rapidly) washed out.

Comparing Figs. 4 and 5 shows that orbital inclinations relax much faster than angular momenta. This is a point already raised in Rauch & Tremaine (1996) (Fig. 2 therein). In particular, they showed in their Section 1.4 that the long-term relaxation of E and L in spherical potential was driven by the NR theory, whereas that of the angular momentum vector L – and thus L_z – was subject to an enhanced relaxation. This out-of-plane relaxation is driven by coherent torques between orbits and is coined vector resonant relaxation (VRR).

4. In-plane diffusion

Assuming that the cluster’s relaxation is driven by local pairwise deflections, its long-term evolution is described by the FP Equation (1). While this equation formally describes the evolution in 3D action space, it is convenient to study relaxation in two-dimensional projections, in particular to improve the signal-to-noise ratio. Yet, such a projection comes at the cost of more intricate theoretical predictions, that require additional integrations along some third action. In this section, following the same approach as in T22, we first focus on in-plane relaxation, that is, relaxation occurring in (J_r, L). Our goal is to compare quantitatively the N-body measurements with the NR prediction.

4.1. N-body measurements

The onset of core collapse is illustrated in Fig. 6 for clusters with N​ = ​10⁵ stars. As in Fig. 3, we recover that rotation does not really impact the rate of core collapse in these linearly stable clusters. Conversely, when the clusters are unstable, such as for (q, α)​ = ​(1, 0.5) (see Fig. 2), we find numerically (not reported in Fig. 6) that the clusters flatten and that relaxation is accelerated.

Fig. 6. Initial evolution of the core radius in clusters with N​ = ​105 stars, as measured in N-body simulations. We use that same convention as in Fig. 2. Interestingly, in these stable clusters, rotation only weakly affects the core contraction.

To probe more precisely relaxation, let us now measure the relaxation rate, ∂F/∂t, in (J_r, L). This is presented in Fig. 7 with details spelled out in Appendix E. Even in the presence of rotation, we recover the same result as in T22, namely all clusters seem to isotropise toward an in-plane DF depending only on energy. Indeed, in radially anisotropic clusters, relaxation depletes radial orbits. And the converse holds for tangentially anisotropic clusters. Importantly, the amount of rotation, α, has no significant impact on the geometry of the in-plane diffusion. Increasing α only (very) weakly accelerates the in-plane relaxation. This is in concordance with the slightly shorter core collapse time observed in Fig. 3.

Fig. 7. Relaxation rate, ∂F/∂t, in (Jr, L) for various velocity anisotropies q (left to right) and rotation parameters α (top and bottom), as measured in N-body simulations. The overall amplitude of the in-plane relaxation rate depends on q, but only very weakly on α.

4.2. Non-resonant prediction

Let us now compute the NR prediction for the relaxation rate in (J_r, L). This prediction was numerically very intensive and required the evaluation of five embedded integrals, as demonstrated in Appendices A and H. Importantly, we stress that in the present case, the LBD (Equation (5)) allowed us to analytically perform the angular part of the orbit average (Equation (4)). This helped reducing the computation time. We refer to Appendix C.1 for details.

The relaxation rate, ∂F/∂t in (J_r, L), presented in Fig. 8 for various anisotropies and rotations, suggests once again an orbital reshuffling toward isotropisation. In addition, increasing the rotation parameter α only very weakly increases the relaxation rate in (J_r, L), while barely impacting the structures observed in the non-rotating case.

Fig. 8. Relaxation rate, ∂F/∂t, in (Jr, L) as in Fig. 7, but here from the NR theory. The prediction matches the N-body measurements from Fig. 7, up to an overall prefactor that depends on q, and (very) weakly on α.

When comparing the N-body measurements (Fig. 7) with the NR prediction (Fig. 8), the NR theory successfully recovers the in-plane relaxation of stable rotating clusters: both figures exhibit strikingly similar structures in action space. We obtain results for the non-rotating case that are in agreement with T22. Yet, even though amplitudes in both approaches are comparable, there is still a (slight) overall prefactor mismatch. The inclusion of rotation has little impact on this mismatch. A more quantitative comparison would require performing (many) more N-body runs, as well as a more precise measurement of the relaxation rate. This would be no light undertaking, and will be the subject of future works.

5. Out-of-plane diffusion

Because they rotate, the present clusters also undergo some out-of-plane diffusion, namely relative to the orbital inclinations cos I. In this section, we focus on relaxation in (J_r, cos I), and we refer to Appendix J for a similar investigation in (L, cos I).

5.1. N-body measurements

The relaxation rates in N-body simulations are presented in Fig. 9. As expected, in the absence of rotation (top row), the relaxation is independent of cos I, and the distribution of orientations remains uniform.

Fig. 9. Relaxation rate, ∂F/∂t, in (Jr, cos I) for various anisotropies q (left to right) and rotation α (from top to bottom), as measured in N-body simulations. The amplitudes and structures observed depend on anisotropy, and show a reshuffling towards isotropisation. Orientations redistribute toward a more affine distribution in cos I. As a result of using binning in action space and finite difference in time to compute ∂F/∂t, the expected initially sharp relaxation occurring at cos I​ = ​0 has been smoothed and is not visible.

Let us now consider the rotating clusters presented in the middle and bottom rows of Fig. 9. In radially anisotropic (q​ = ​1) and isotropic (q​ = ​0) clusters, the systems lose stars in the prograde region (cos I​ > ​0) and gain stars in the retrograde region (cos I​ < ​0). As one increases α (i.e. as one increases the net rotation), this trend strengthens. This is in agreement with the orbital reshuffling observed in Fig. 8, where the clusters always tend to isotropise their in-plane distribution. Finally, we note that the highest diffusion rates in inclination are observed near cos I​ = ​0. This corresponds precisely to the location of the discontinuity of the LBD.

At first glance, the tangentially anisotropic case (q​ = ​ − 6), as given by the right column in Fig. 9, may seem different. Indeed, in that case a depletion of orbits for small J_r is observed, whatever cos I. Conversely, the number of orbits systematically increases for large J_r, whatever cos I. This is directly linked to the in-plane isotropisation of the cluster (Fig. 8). Indeed, these clusters being tangentially biased, the in-plane diffusion occurs towards higher J_r. Overall, all clusters, independently of their anisotropies, evolve towards smoother distribution of orbital inclinations.

As a complement of Fig. 9, Appendix J also presents the relaxation rates in (L, cos I). It reaches similar results as for diffusion in (J_r, cos I). Namely, the in-plane distributions tend to isotropise and orbital inclinations diffuse so as to reach smoother distributions.

In addition, we stress that due to the use of binning and finite differentiation to compute ∂F/∂t, the expected sharp evolution occurring at cos I​ = ​0 has been smoothed out. Furthermore, the amplitude of relaxation depends on the time interval used for finite differentiation and on the bin size. While this does not drastically change the observations away from cos I​ = ​0, any measurement in the neighborhood of cos I​ = ​0 is tricky to perform: it should be taken with appropriate caution.

5.2. Non-resonant prediction

Let us now perform the same investigation in (J_r, cos I) using the NR theory. Obtaining a satisfying prediction near cos I​ = ​0 required finely sampling the NR orbital integrals, as detailed in Appendix H.2. The NR prediction is illustrated in Fig. 10.

Fig. 10. Relaxation rate, ∂F/∂t, in (Jr, cos I) as in Fig. 9, but here from the NR theory (see Appendix H.2 for details). While the in-plane isotropisation and the smoothing of the inclination distribution is recovered, finer details disagree. Here, the amount of relaxation was estimated using finite differentiation using separately regions with positive or negative cos I In effect, this removes the δ′(cos I) singularity in cos I​ = ​0, as discussed in Appendix I. Importantly, the removal of this singularity does not affect the prediction away from cos I = 0. In addition, because this singularity is diluted in N-body measurements – as discussed in Fig. 9 – we do not represent the NR prediction here.

First, in the absence of rotation, the NR theory predicts a relaxation rate independent of cos I: this is a reassuring sanity check. Yet, when comparing the top rows of Figs. 9 and 10, we note some discrepancies. Namely, the locations of the line of zero relaxation (dashed lines) slightly differ.

The difference is much more striking in the presence of rotation (middle and bottom rows of Fig. 10). In both figures, the orbital reshuffling still operates, and the clusters redistribute their orbits towards smoother distributions of inclinations. Yet, the diffusion structures predicted by the NR theory do not align with the N-body ones. Indeed, N-body simulations (Fig. 9) present “round” structures as one goes from cos I​ = ​0 to cos I​ = ​1. In addition, the relaxation rate decreases as one considers orbits within to rotation plane, that is, as one considers cos I​ → ​ ± 1. This is in sharp contrast with the “straighter” structures predicted by the NR theory (Fig. 10). Nevertheless, one should keep in mind that Fig. 10 cuts out the singularity at cos I​ = ​0. The comparison to the N-body measurements near cos I​ = ​0 is made difficult by the fact that diffusion at early time is very fast in that region, as shown by Fig. K.4. As a consequence, the measured relaxation rate may suffer from truncation errors – due to, for example, the finite differentiation scheme – the closer one gets to cos I​ = ​0. However, as Fig. K.4 shows, the initial discontinuity appears to have little impact on diffusion both far away from cos I​ = ​0 and beyond t​ = ​0.1 t_rh. The obtention of a more robust measure would require both a finer time difference and a finer binning of action space. These two improvements would require an ensemble-average over much more realisations to reduce dispersion. A similar conclusion is reached in Appendix J when considering relaxation in (L, cos I).

Before concluding, let us elaborate further on the failure of the NR theory to model out-of-plane relaxation. For α​ ≠ ​0, recall that cos I​ = ​ ± 1 corresponds to orbits within the rotation plane, while cos I​ = ​0 corresponds to orbits perpendicular to the rotation plane. A test star orbiting in this rotating cluster will torque with all the other stars of the cluster. Equivalently, this means that this star will typically be subject to a torque with the total angular momentum of the rotating cluster, ⟨L⟩​ ∝ ​α e_z, the strength of which depends on the test star’s orbital orientation. This is what we observe in Fig. 9. However, the NR theory only takes into account local deflections, and therefore cannot take these coherent interactions into account.

Following the work of Meiron & Kocsis (2019) (see also references therein), we anticipate for VRR to play a key role in these coherent interactions – even after factoring the truncation errors due to the finite difference scheme and the singularity at cos I​ = ​0, both discussed earlier in this section. Indeed, VRR, which is driven by persistent torques between orbital planes, might be needed to explain (i) the faster out-of-plane relaxation (Section 3.3); (ii) the discrepancies between the N-body measurements and the NR predictions. Investigating quantitatively the efficiency of VRR in these systems will be the topic of future investigations.

6. Conclusions and perspectives

6.1. Conclusion

In this paper, we used the NR formalism to probe the long-term time-evolution of rotating anisotropic globular clusters. First, using N-body simulations, we showed how rotation only has a weak impact on the time of core collapse. As such, in contrast with prior research (see, e.g., Hachisu 1979), we did not observe any gravo-gyro catastrophe that could expedite core collapse (see, e.g., Einsel & Spurzem 1999). Since we introduced rotation in a different fashion, this does not contradict previous results.

Focusing on in-plane diffusion, we showed how the NR prediction successfully recovers all the intricacies of the N-body measurements. This is in line with the non-rotating results from T22. Yet, although the diffusion structures in action space closely align, there is still an overall amplitude mismatch between the N-body measurements and the kinetic prediction.

We subsequently turned our interest on out-of-plane relaxation, that is, the redistribution of orbital inclinations. We pointed out the similarities and differences between the N-body measurements and the NR prediction. In both approaches, we observe a systematic reshuffling of orbital inclinations from overpopulated regions to underpopulated ones. As such, the distribution of inclinations get smoother through time, and this process is accelerated by rotation. On the one hand, we observed that the initial discontinuity at cos I​ = ​0 induced a fast relaxation at small times. Surely, this strongly affects any measurement of the relaxation rate in that region. On the other hand, for orbits aligned with the cluster’s rotation plane – where diffusion is much slower and measurements are much easier – the local NR theory does not match with the N-body measurements. We argue that this mismatch originates from the fact that the NR formalism is sourced by local deflections, which neglect the coherent torques between the orbital planes. Accounting for the (efficient) contribution of persistent torques requires the use of VRR.

6.2. Perspectives

Here, the NR formalism was used to describe the very onset of relaxation. Yet, the short duration of our simulations and our limited sample of initial conditions prevented us from fully assessing the asymptotic late time orbital distribution. Naturally, it would be interesting to push the N-body integrations further in time. While some studies carried out such long-term N-body simulations (Tiongco et al. 2020; Livernois et al. 2022), orbital inclinations were, unfortunately, not their focus. Along the same line, it would also be of interest to integrate the 3D FP equation itself. In isotropic clusters, this process has been conducted, as demonstrated for instance in the work of Vasiliev (2015). However, the scenario involving anisotropic, rotating systems has yet to be comprehensively investigated.

Ultimately, if rotation is sufficiently large, it may induce a flattening of the cluster. In that case, Stäckel systems (Dejonghe & de Zeeuw 1988) could be used to model realistic flattened rotating structures. In particular, this would still ensure the integrability of the mean potential, that is the existence of explicit angle-action coordinates. Tailoring secular theory to such setups will be the topic of future works.

In Section 5, we pointed out how the NR theory does not reproduce the out-of-plane diffusion observed in N-body simulations. To understand this mismatch, it would be valuable to examine the predictions of the resonant relaxation (RR) formalism (see, e.g., Hamilton et al. 2018; Fouvry et al. 2021) in the presence of rotation. While the inhomogeneous Landau equation was already investigated in non-rotating isotropic clusters (Hamilton et al. 2018; Fouvry et al. 2021), its implementation in rotating spheres with an explicit dependence on L_z should be the next step. Furthermore, to account for collective effects, an explicit implementation of the Balescu–Lenard equation (Heyvaerts 2010) might prove necessary to match the details of stacked N-body measurements.

Similarly, to get a better handle on the mismatch reported in section 5, one should investigate quantitatively the efficiency of VRR in these clusters (Rauch & Tremaine 1996). This should allow us to determine if indeed coherent and long-lasting torques between orbital planes are the missing ingredients to describe the relaxation of orientations. While this has been extensively studied in galactic nuclei (see, e.g., Eilon et al. 2009; Kocsis & Tremaine 2011, 2015; Szölgyén & Kocsis 2018; Fouvry et al. 2019; Szölgyén et al. 2019, 2021; Magnan et al. 2022), the implementation of VRR in globular clusters should be the topic of future work (see Meiron & Kocsis 2019, for a preliminary investigation).

It would be interesting to leverage multiple masses into the theory, so as to capture (radial or inclination) segregation effects (Meiron & Kocsis 2019). Dekel et al. (2023) argues for instance that the excess of massive galaxies found at very high redshift by the James Webb Space Telescope could originate from early feedback-free star formation within dense, possibly rotating stellar clusters made of massive stars. The fate of such clusters should be captured by the present secular theory. In particular, it could lead to the formation of intermediate mass BHs (Greene & Ho 2004; Greene et al. 2020) and possibly seeds for supermassive BHs (Kormendy & Ho 2013).

Data availability

The data underlying this article is available through reasonable request to the authors. The code, written in JULIA (Bezanson et al. 2017), computing the NR diffusion coefficients in anisotropic rotating clusters is available at the URL: https://github.com/KerwannTEP/CARP.

Acknowledgments

This work is partially supported by the grant https://www.secular-evolution.org ANR-19-CE31-0017 of the French Agence Nationale de la Recherche, and by the Idex Sorbonne Université. We also thank the KITP for hosting the workshop https://www.cosmicweb23.org during which this project was advanced. This work is partially supported by the National Science Foundation under Grant No. NSF PHY-1748958 and Grant No. AST-2310362. We are grateful to M. Roule, M. Petersen, A. L. Varri and the members of the Segal collaboration for numerous suggestions during the completion of this work. We thank Stéphane Rouberol for the smooth running of the Infinity cluster, where the simulations were performed. Finally, we thank our referee, Prof. Douglas Heggie, for his remarks, which contributed greatly to the improvement of this paper.

References

Appendix A: Non-resonant theory

In this Appendix, we detail the derivation of equation (3). This is the backbone of our computation of the NR theory.

A.1. Local velocity coefficients

The starting point are equations (2). We consider this relation within the frame from Fig. A.1. The local velocity deflections are then given by

$$\begin{array}{cc}\hfill ⟨\mathrm{\Delta }{v}_{\Vert }⟩& =⟨\mathrm{\Delta }{v}_1⟩,\hfill \end{array}$$(A.1a)

$$\begin{array}{cc}\hfill ⟨{(\mathrm{\Delta }{v}_{\Vert })}^2⟩& =⟨{(\mathrm{\Delta }{v}_1)}^2⟩,\hfill \end{array}$$(A.1b)

$$\begin{array}{cc}\hfill ⟨{(\mathrm{\Delta }{v}_{\perp })}^2⟩& =⟨{(\mathrm{\Delta }{v}_2)}^2⟩+⟨{(\mathrm{\Delta }{v}_3)}^2⟩,\hfill \end{array}$$(A.1c)

Fig. A.1. Tailored frame used to compute the parallel and perpendicular local velocity deflections in equations (A.1). By construction, the test star’s angular momentum L is along the axis e2.

where the coordinate 1 is the z coordinate, parallel to v. The coordinates 2 and 3 are defined such that the projection of r on the (Oxy) plane is along e₃ (Fig. A.1). We now introduce w​ = ​v​ − ​v′ and change the integration variables from v′ to w in equation (2). We get

$$\begin{array}{cc}\hfill ⟨\mathrm{\Delta }{v}_i⟩& =-2A{\displaystyle \int \mathrm{d}}\mathit{w}\frac{w_i}{w^3}{F}_{\mathrm{rot}},\hfill \end{array}$$(A.2a)

$$\begin{array}{cc}\hfill ⟨\mathrm{\Delta }{v}_i\mathrm{\Delta }{v}_j⟩& =A{\displaystyle \int \mathrm{d}}\mathit{w}\frac{w^2{\delta }_{\mathit{ij}}-w_i{w}_j}{w^3}{F}_{\mathrm{rot}},\hfill \end{array}$$(A.2b)

with the shortened notation F_rot​ = ​F_rot(r, v′​ = ​v​ − ​w). Using these equations allows us to compute the needed coefficients from equation (A.1) through

$$\begin{array}{c}\hfill \left[\begin{array}{c}⟨\mathrm{\Delta }{v}_1⟩\\ ⟨{(\mathrm{\Delta }{v}_1)}^2⟩\\ ⟨{(\mathrm{\Delta }{v}_2)}^2⟩\\ ⟨{(\mathrm{\Delta }{v}_3)}^2⟩\end{array}\right]=A{\displaystyle \int \frac{\mathrm{d}\mathit{w}}{w^3}}\left[\begin{array}{c}-2w_1\\ w^2-w_1^2\\ w^2-w_2^2\\ w^2-w_3^2\end{array}\right]F_{\mathrm{rot}}.\end{array}$$(A.3)

We now define the spherical coordinates associated with this frame (Fig. A.1)

$$\begin{array}{c}\hfill w_1=wcos\vartheta ;w_2=wsin\vartheta cos\varphi ;w_3=wsin\vartheta sin\varphi .\end{array}$$(A.4)

Injecting equations (A.4) into equations (A.3), we obtain

$$\begin{array}{c}\hfill \left[\begin{array}{c}⟨\mathrm{\Delta }{v}_1⟩\\ ⟨{(\mathrm{\Delta }{v}_1)}^2⟩\\ ⟨{(\mathrm{\Delta }{v}_2)}^2⟩\\ ⟨{(\mathrm{\Delta }{v}_3)}^2⟩\end{array}\right]=A{\displaystyle \int \mathrm{d}}w\mathrm{d}\vartheta \mathrm{d}\varphi sin\vartheta \left[\begin{array}{c}-2cos\vartheta \\ w{sin}^2\vartheta \\ w(1-{sin}^2\vartheta {cos}^2\varphi )\\ w(1-{sin}^2\vartheta {sin}^2\varphi )\end{array}\right]F_{\mathrm{rot}}.\end{array}$$(A.5)

Here, the background distribution, F_rot, is to be evaluated in $F_{\text{rot}}(\mathit{r},{\mathit{v}}^{\prime })\text{}=\text{}{F}_{\text{rot}}(E^{\prime },L^{\prime },L_z^{\prime })$. We compute these arguments in appendix A.2. Finally injecting equation (A.5) into equation (A.1), we obtain our main result, namely equation (3), as given in the main text.

A.2. Arguments of the background distribution

Let us now compute the arguments E′, L′ and $L_z^{\prime }$ for the background distribution function, $F_{\text{rot}}(E^{\prime },L^{\prime },L_z^{\prime })$, in equation (2).

We start with E′. It reads

$$\begin{array}{cc}\hfill E^{\prime }& =\psi (r)+\frac{{\mathit{v}}^{\prime 2}}{2}=\psi (r)+\frac{{(\mathit{v}-\mathit{w})}^2}{2}=E+\frac{w^2}{2}-wvcos\vartheta ,\hfill \end{array}$$(A.6)

since v​ = ​v e₁ (see Fig. A.1). This form allows us to get an upper bound on the w-integral of equation (A.5), above which E′​ > ​0. As soon as E′​ > ​0, the background star is unbound and its DF vanishes.

The computation of L′ and $L_z^{\prime }$ is best achieved by using a different coordinate system, as illustrated in Figs. A.1 and A.2. More precisely, the frame (1, 2, 3) is related to the (X, Y, Z) one through the relations

$$\begin{array}{c}\hfill \left[\begin{array}{c}{\mathbf{e}}_1\hfill \\ {\mathbf{e}}_2\hfill \\ {\mathbf{e}}_3\hfill \end{array}\right]=\left[\begin{array}{ccc}\frac{v_{\mathrm{r}}}{v}cos\phi -\frac{v_{\mathrm{t}}}{v}sin\phi & \frac{v_{\mathrm{r}}}{v}sin\phi +\frac{v_{\mathrm{t}}}{v}cos\phi & 0\\ 0& 0& 1\\ \frac{v_{\mathrm{r}}}{v}sin\phi +\frac{v_{\mathrm{t}}}{v}cos\phi & \frac{v_{\mathrm{t}}}{v}sin\phi -\frac{v_{\mathrm{r}}}{v}cos\phi & 0\end{array}\right]\left[\begin{array}{c}{\displaystyle {\mathbf{e}}_{\mathbf{X}}}\hfill \\ {\displaystyle {\mathbf{e}}_{\mathbf{Y}}}\hfill \\ {\displaystyle {\mathbf{e}}_Z}\hfill \end{array}\right],\end{array}$$(A.7)

Fig. A.2. Illustration of the frames used to compute the argument of the background DF, $F_{\text{rot}}(E^{\prime },L^{\prime },L_z^{\prime })$. In the left panel, the cluster is seen from the side of the rotation plane, (Oxy). In the right panel, the (OXY)-plane is the orbital plane of the test star, seen from above.

with v_r (resp. v_t) the radial (resp. tangential) component of the test star’s velocity. We can finally relate the frame (X, Y, Z) to the frame (x, y, z) through the relations

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\mathbf{e}}_X\hfill \\ {\mathbf{e}}_Y\hfill \\ {\mathbf{e}}_Z\hfill \end{array}\right]& =\left[\begin{array}{ccc}cosI& 0& sinI\\ 0& 1& 0\\ -sinI& 0& cosI\end{array}\right]\left[\begin{array}{c}{\displaystyle {\mathbf{e}}_{\mathbf{x}}}\hfill \\ {\displaystyle {\mathbf{e}}_{\mathbf{y}}}\hfill \\ {\displaystyle {\mathbf{e}}_{\mathbf{z}}}\hfill \end{array}\right].\hfill \end{array}$$(A.8)

We can now compute L, the norm of L. As it is a norm, we can choose any frame to compute it. The frame (1, 2, 3) is a convenient choice as this frame is used to compute the integrands of equations (A.3). We write

$$\begin{array}{cc}\hfill L^{\prime 2}& ={(r_1{v}_2^{\prime }-r_2{v}_1^{\prime })}^2+{(r_1{v}_3^{\prime }-r_3{v}_1^{\prime })}^2+{(r_2{v}_3^{\prime }-r_3{v}_2^{\prime })}^2,\hfill \end{array}$$(A.9)

where (r₁, r₂, r₃) are the coordinates of r in the frame (1, 2, 3). Using equations (A.7) and (A.8), we obtain

$$\begin{array}{c}\hfill r_1=\frac{v_{\mathrm{r}}}{v}r;r_2=0;r_3=\frac{v_{\mathrm{t}}}{v}r.\end{array}$$(A.10)

Therefore, one gets

$$\begin{array}{cc}\hfill L^{\prime }& =r\sqrt{v_2^{\prime 2}+(\frac{v_{\mathrm{r}}}{v}{v}_3^{\prime }-\frac{v_{\mathrm{t}}}{v}{v}_1^{\prime }{)}^2}\hfill \\ \hfill & =r\sqrt{{(wsin\vartheta cos\varphi )}^2+(v_{\mathrm{t}}+\frac{v_{\mathrm{r}}}{v}wsin\vartheta sin\varphi -\frac{v_{\mathrm{t}}}{v}wcos\vartheta {)}^2},\hfill \end{array}$$(A.11)

given that v′​ = ​v − w.

Let us finally compute $L_z^{\prime }$. In the (x, y, z) frame, it reads

$$\begin{array}{c}\hfill L_z^{\prime }=x{v}_y^{\prime }-y{v}_x^{\prime }.\end{array}$$(A.12)

From equation (A.8), the positions read

$$\begin{array}{cc}\hfill x& =XcosI=rcos\phi cosI,\hfill \end{array}$$(A.13a)

$$\begin{array}{cc}\hfill y& =Y=rsin\phi ,\hfill \end{array}$$(A.13b)

$$\begin{array}{cc}\hfill z& =XsinI=rcos\phi sinI,\hfill \end{array}$$(A.13c)

using X​ = ​r cos φ and Y​ = ​r sin φ. We then obtain the velocities

$$\begin{array}{cc}\hfill v_x& =v_{\mathrm{r}}cos\phi cosI-v_{\mathrm{t}}sin\phi cosI,\hfill \end{array}$$(A.14a)

$$\begin{array}{cc}\hfill v_y& =v_{\mathrm{r}}sin\phi +v_{\mathrm{t}}cos\phi ,\hfill \end{array}$$(A.14b)

$$\begin{array}{cc}\hfill v_z& =v_{\mathrm{r}}cos\phi sinI-v_{\mathrm{t}}sin\phi sinI,\hfill \end{array}$$(A.14c)

by taking the time derivative of equations (A.13). Let us now compute the background velocity components. First, using the left panel of Fig. A.2, we have

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{v}_x^{\prime }\hfill \\ v_y^{\prime }\hfill \\ v_z^{\prime }\hfill \end{array}\right]& =\left[\begin{array}{ccc}cosI& 0& -sinI\\ 0& 1& 0\\ sinI& 0& cosI\end{array}\right]\left[\begin{array}{c}{\displaystyle v_X^{\prime }}\hfill \\ {\displaystyle v_Y^{\prime }}\hfill \\ {\displaystyle v_Z^{\prime }}\hfill \end{array}\right].\hfill \end{array}$$(A.15)

Second, using the relation in equations (A.7) (Fig. A.1 and right panel of Fig. A.2), we have

$$\begin{array}{c}\hfill \left[\begin{array}{c}{v}_X^{\prime }\hfill \\ v_Y^{\prime }\hfill \\ v_Z^{\prime }\hfill \end{array}\right]=\left[\begin{array}{ccc}\frac{v_{\mathrm{r}}}{v}cos\phi -\frac{v_{\mathrm{t}}}{v}sin\phi & 0& \frac{v_{\mathrm{r}}}{v}sin\phi +\frac{v_{\mathrm{t}}}{v}cos\phi \\ \frac{v_{\mathrm{r}}}{v}sin\phi +\frac{v_{\mathrm{t}}}{v}cos\phi & 0& \frac{v_{\mathrm{t}}}{v}sin\phi -\frac{v_{\mathrm{r}}}{v}cos\phi \\ 0& 1& 0\end{array}\right]\left[\begin{array}{c}{\displaystyle v_1^{\prime }}\hfill \\ {\displaystyle v_2^{\prime }}\hfill \\ {\displaystyle v_3^{\prime }}\hfill \end{array}\right].\end{array}$$(A.16)

Overall, we can rewrite equation (A.12) as

$$\begin{array}{cc}\hfill L_z^{\prime }& =r(\frac{v_1^{\prime }{v}_{\mathrm{t}}}{v}-\frac{v_3^{\prime }{v}_{\mathrm{r}}}{v})cosI+r{v}_2^{\prime }sin\phi sinI\hfill \\ \hfill & =r(v_{\mathrm{t}}-\frac{v_{\mathrm{t}}}{v}wcos\vartheta +\frac{v_{\mathrm{r}}}{v}wsin\vartheta sin\varphi )cosI\hfill \\ \hfill & -rwsin\vartheta cos\varphi sin\phi sinI.\hfill \end{array}$$(A.17)

Finally, using the Cauchy–Schwarz inequality, one can check that

$$\begin{array}{c}\hfill |L_z^{\prime }|\le L^{\prime }\sqrt{{cos}^{2}I+{sin}^2\phi {sin}^{2}I},\end{array}$$(A.18)

which ensures that $\left|L_z^{\prime }\right|\text{}\le \text{}{L}^{\prime }$.

Equations (A.6), (A.11) and (A.17) are the main results of this section. They provide us with explicit expressions for the arguments at which to evaluate $F_{\text{rot}}\text{}=\text{}{F}_{\text{rot}}(E^{\prime },L^{\prime },L_z^{\prime })$ in equation (3).

Appendix B: 3D diffusion coefficients

The local in-plane diffusion coefficients, that is, the diffusion coefficients in E and L, are given in appendix C of Bar-Or & Alexander (2016). We reproduce them here for completeness. They read

$$\begin{array}{cc}\hfill ⟨\mathrm{\Delta }E⟩& =\frac{1}{2}⟨{(\mathrm{\Delta }{v}_{\Vert })}^2⟩+\frac{1}{2}⟨{(\mathrm{\Delta }{v}_{\perp })}^2⟩+v⟨\mathrm{\Delta }{v}_{\Vert }⟩,\hfill \end{array}$$(B.1a)

$$\begin{array}{cc}\hfill ⟨{(\mathrm{\Delta }E)}^2⟩& =v^2⟨{(\mathrm{\Delta }{v}_{\Vert })}^2⟩,\hfill \end{array}$$(B.1b)

$$\begin{array}{cc}\hfill ⟨\mathrm{\Delta }L⟩& =r\frac{v_{\mathrm{t}}}{v}⟨\mathrm{\Delta }{v}_{\Vert }⟩+\frac{r^2}{4L}⟨{(\mathrm{\Delta }{v}_{\perp })}^2⟩,\hfill \end{array}$$(B.1c)

$$\begin{array}{cc}\hfill ⟨(\mathrm{\Delta }{L}^2)⟩& =r^2\frac{v_{\mathrm{t}}^2}{v^2}⟨{(\mathrm{\Delta }{v}_{\Vert })}^2⟩+\frac{r^2}{2}\frac{v_{\mathrm{r}}^2}{v^2}⟨{(\mathrm{\Delta }{v}_{\perp })}^2⟩,\hfill \end{array}$$(B.1d)

$$\begin{array}{cc}\hfill ⟨\mathrm{\Delta }E\mathrm{\Delta }L⟩& =L⟨{(\mathrm{\Delta }{v}_{\Vert })}^2⟩.\hfill \end{array}$$(B.1e)

In order to compute the 3D FP equation (1), we also need the local diffusion coefficients involving L_z​ = ​L​ ⋅ ​e_z. Performing a first-order variation of L​ = ​r × v, we can write

$$\begin{array}{cc}\hfill \mathrm{\Delta }\mathit{L}& =(\frac{\mathrm{\Delta }{v}_{\Vert }}{v}-\frac{v_{\mathrm{r}}}{v_{\mathrm{t}}}\frac{\mathrm{\Delta }{v}_3}{v})\mathit{L}+L\frac{\mathrm{\Delta }{v}_2}{v}(\frac{v_{\mathrm{r}}}{v_{\mathrm{t}}}{\mathbf{e}}_3-\widehat{\mathit{v}}),\hfill \end{array}$$(B.2a)

$$\begin{array}{cc}\hfill {\mathbf{e}}_3& =\frac{\mathit{v}\times \mathit{L}}{|\mathit{v}\times \mathit{L}|}=\frac{v\widehat{\mathit{r}}-v_{\mathrm{r}}\widehat{\mathit{v}}}{v_{\mathrm{t}}},\hfill \end{array}$$(B.2b)

where (e₁, e₂, e₃) is the frame illustrated in Fig. A.1. We can then write

$$\begin{array}{c}\hfill {\mathbf{e}}_3·{\mathbf{e}}_z=sinI(\frac{v_{\mathrm{r}}}{v}sin\phi +\frac{v_{\mathrm{t}}}{v}cos\phi ).\end{array}$$(B.3)

As a result, it follows that

$$\begin{array}{c}\hfill \frac{v_{\mathrm{r}}}{v}{\mathbf{e}}_3·{\mathbf{e}}_z-\frac{v_z{v}_{\mathrm{t}}}{v^2}=sin\phi sinI,\end{array}$$(B.4)

with φ given in Fig. (A.2) and v_z by equation (A.14c). Overall, we finally have

$$\begin{array}{cc}\hfill ⟨\mathrm{\Delta }{L}_z⟩& =\frac{L_z}{v}⟨\mathrm{\Delta }{v}_{\Vert }⟩,\hfill \end{array}$$(B.5a)

$$\begin{array}{cc}\hfill ⟨(\mathrm{\Delta }{L}_z^2⟩& =(\frac{L_z}{L}{)}^2(\frac{L^2}{v^2}⟨{(\mathrm{\Delta }{v}_{\Vert })}^2⟩+\frac{1}{2}\frac{r^2{v}_{\mathrm{r}}^2}{v^2}⟨{(\mathrm{\Delta }{v}_{\perp })}^2⟩)\hfill \end{array}$$(B.5b)

$$\begin{array}{cc}& +\frac{r^2{sin}^2\phi }{2}(1-\frac{L_z^2}{L^2})⟨{(\mathrm{\Delta }{v}_{\perp })}^2⟩,\hfill \\ \hfill ⟨\mathrm{\Delta }E\mathrm{\Delta }{L}_z⟩& =L_z⟨{(\mathrm{\Delta }{v}_{\Vert })}^2⟩,\hfill \end{array}$$(B.5c)

$$\begin{array}{cc}\hfill ⟨\mathrm{\Delta }L\mathrm{\Delta }{L}_z⟩& =\frac{L_z}{L}(\frac{L^2}{v^2}⟨{(\mathrm{\Delta }{v}_{\Vert })}^2⟩+\frac{1}{2}\frac{r^2{v}_{\mathrm{r}}^2}{v^2}⟨{(\mathrm{\Delta }{v}_{\perp })}^2⟩).\hfill \end{array}$$(B.5d)

Equations (B.1) and (B.5) are the main results of this section. They allow us to compute the local diffusion coefficients in (E, L, L_z) from the local velocity diffusion coefficients given by equation (3).

Appendix C: Orbit average

The computation of the global diffusion coefficients in the FP equation (1) requires that we orbit average the local diffusion coefficients over the mean field orbit of the test star. Since the background cluster is rotating, this average involves a radial integration, but also an angular one. This is what we detail here.

C.1. Angular orbit average

We start from the generic formula of equation (4). There, the average with respect to φ must be performed with care since $L_z^{\prime }$ depends on φ (equation A.17). Fortunately, taking advantage of the particular structure of the diffusion coefficients in E, L and L_z (see equations B.1 and B.5), we only need to compute the following four integrals over φ

$$\begin{array}{cc}& {\displaystyle \int \frac{\mathrm{d}\phi }{2\pi }}⟨\mathrm{\Delta }{v}_{\Vert }⟩,{\displaystyle \int \frac{\mathrm{d}\phi }{2\pi }}⟨{(\mathrm{\Delta }{v}_{\Vert })}^2⟩,\hfill \\ \hfill & {\displaystyle \int \frac{\mathrm{d}\phi }{2\pi }}⟨{(\mathrm{\Delta }{v}_{\perp })}^2⟩,{\displaystyle \int \frac{\mathrm{d}\phi }{2\pi }}cos2\phi ⟨{(\mathrm{\Delta }{v}_{\perp })}^2⟩,\hfill \end{array}$$(C.1)

where we used ${sin}^2\phi =\frac{1}{2}(1-cos2\phi )$ for convenience.

In equation (3), the only dependence with respect to φ is in the dependence with respect to $L_z^{\prime }$ in F_rot. In the particular case of the LBD (equation 5), we can perform explicitly the average over φ. Computing the integrals (C.1) only requires the evaluation of the two non-trivial integrals

$$\begin{array}{cc}\hfill A_1& ={\displaystyle {\int }_0^{2\pi }}\frac{\mathrm{d}\phi }{2\pi }\mathrm{sgn}(L_z^{\prime }/L^{\prime }),\hfill \end{array}$$(C.2a)

$$\begin{array}{cc}\hfill A_2& ={\displaystyle {\int }_0^{2\pi }}\frac{\mathrm{d}\phi }{2\pi }cos2\phi \mathrm{sgn}(L_z^{\prime }/L^{\prime }),\hfill \end{array}$$(C.2b)

with $L_z^{\prime }$ following from equation (A.17). Introducing $\zeta \text{}=\text{}-(v_1^{\prime }{v}_{\text{t}}/v\text{}-\text{}{v}_3^{\prime }{v}_{\text{r}}/v)\cos I$ and $\mu \text{}=\text{}\zeta /(v_2^{\prime }\sin I)$, these two integrals can be explicitly computed. Their expressions are gathered in Table C.1. This is the main result of this section.

From the numerical point of view, this analytical calculation allows us to speed up the computation of the NR predictions, as the angular integral is removed. In addition, it also explicitly deals with the discontinuity of the LBD DF: this enhances the numerical stability.

C.2. Radial orbit average

Having performed the angle average in Appendix C.1, equation (4) takes the form

$$\begin{array}{c}\hfill D_X=\frac{{\mathrm{\Omega }}_r}{\pi }{\displaystyle {\int }_{r_{\mathrm{p}}}^{r_{\mathrm{a}}}}\frac{\mathrm{d}r}{|v_{\mathrm{r}}|}{⟨\mathrm{\Delta }X⟩}_{\phi }(r),\end{array}$$(C.3)

where ⟨ ⋅ ⟩_φ stands for the angle average. As described in appendix F2 of T22, equation (C.3) can be efficiently performed by changing variable from the radius to an effective anomaly. In particular, this removes the integrable singularity 1/|v_r| at pericentre and apocentre. We follow the exact same approach here, and perform the radial average using 50 sampling nodes.

Appendix D: Fokker–Planck and coordinates

Depending on the quantities we wish to investigate, we may wish to change of coordinates system. Fortunately, it is possible to transform one FP description into another under a change of variable. We describe a few relevant examples in this appendix.

D.1. Generic change of coordinates

Given some coordinates x, the FP equation generically reads (Risken 1996)

$$\begin{array}{cc}\hfill \frac{\partial F(\mathit{x},t)}{\partial t}& =-\frac{\partial }{\partial \mathit{x}}·{\mathcal{F}}_{\mathit{x}}(\mathit{x})\hfill \\ \hfill & =-\frac{\partial }{\partial \mathit{x}}·[{\mathit{D}}_{\mathit{x}}(\mathit{x})F(\mathit{x})-\frac{1}{2}\frac{\partial }{\partial \mathit{x}}·({\mathit{D}}_{\mathit{x}\mathit{x}}(\mathit{x})F(\mathit{x}))],\hfill \end{array}$$(D.1)

with the diffusion coefficients

$$\begin{array}{c}\hfill {\mathit{D}}_{\mathit{x}}(\mathit{x})=\left[\begin{array}{c}{D}_{x_1}\\ D_{x_2}\\ D_{x_3}\end{array}\right];{\mathit{D}}_{\mathit{x}\mathit{x}}(\mathit{x})=\left[\begin{array}{ccc}{D}_{x_1{x}_1}& D_{x_1{x}_2}& D_{x_1{x}_3}\\ D_{x_1{x}_2}& D_{x_2{x}_2}& D_{x_2{x}_3}\\ D_{x_1{x}_3}& D_{x_2{x}_3}& D_{x_3{x}_3}\end{array}\right],\end{array}$$(D.2)

and the DF in x space, F(x). In our case, depending on the context, x​ = ​(x₁, x₂, x₃) may either stand for (E, L, L_z), (J_r, L, L_z) or (J_r, L, cos I).

Following Risken (1996) and Bar-Or & Alexander (2016), one can easily rewrite the FP equation within some new coordinates x′(x). The new diffusion coefficients read

$$\begin{array}{cc}\hfill D_{x_l^{\prime }}& ={\displaystyle \sum _k}\frac{\partial x_l^{\prime }}{\partial x_k}{D}_{x_k}+\frac{1}{2}{\displaystyle \sum _{k,r}}\frac{{\partial }^2{x}_l^{\prime }}{\partial x_r\partial x_k}{D}_{x_r{x}_k},\hfill \end{array}$$(D.3a)

$$\begin{array}{cc}\hfill D_{x_l^{\prime }{x}_m^{\prime }}& ={\displaystyle \sum _{k,r}}\frac{\partial x_l^{\prime }}{\partial x_r}\frac{\partial x_m^{\prime }}{\partial x_k}{D}_{x_r{x}_k}.\hfill \end{array}$$(D.3b)

These coefficients source a FP equation in x′-space for the DF, F′(x′), reading

$$\begin{array}{c}\hfill F^{\prime }({\mathit{x}}^{\prime })=|\frac{\partial \mathit{x}}{\partial {\mathit{x}}^{\prime }}|F(\mathit{x}),\end{array}$$(D.4)

with |∂x/∂x′| the inverse Jacobian of the coordinate transform.

D.2. From (E, L, L_z) to (J_r, L, L_z)

Following section 2, we have at our disposal diffusion coefficients within the coordinates (E, L, L_z). Owing to equations (D.3), the diffusion coefficients in the action space (J_r[E, L],L, L_z) are easily computed. This is already detailed in appendix F3 of T22. Ultimately, we have at our disposal the diffusion coefficients

$$\begin{array}{c}\hfill {\mathit{D}}_1(\mathit{J})=\left[\begin{array}{c}{D}_{J_{\mathrm{r}}}\hfill \\ D_L\hfill \\ D_{L_z}\hfill \end{array}\right];{\mathit{D}}_2(\mathit{J})=\left[\begin{array}{ccc}{D}_{J_{\mathrm{r}}{J}_{\mathrm{r}}}& D_{J_{\mathrm{r}}L}& D_{J_{\mathrm{r}}{L}_z}\\ D_{J_{\mathrm{r}}L}& D_{\mathit{LL}}& D_{L{L}_z}\\ D_{J_{\mathrm{r}}{L}_z}& D_{L{L}_z}& D_{L_z{L}_z}\end{array}\right],\end{array}$$(D.5)

which source the FP evolution of F_rot(J_r, L, L_z), the cluster’s DF in (J_r, L, L_z).

D.3. From (J_r, L, L_z) to (J_r, L, cos I)

Let us now change of coordinates from J​ = ​(J_r, L, L_z) to J_c​ = ​(J_r, L, cos I​ = ​L_z/L). Applying equations (D.3), we obtain the new diffusion coefficients

$$\begin{array}{cc}\hfill D_{cosI}& =-\frac{cosI{D}_L}{L}+\frac{D_{L_z}}{L}+\frac{cosI{D}_{\mathit{LL}}}{L^2}-\frac{D_{L{L}_z}}{L^2},\hfill \end{array}$$(D.6a)

$$\begin{array}{cc}\hfill D_{J_{\mathrm{r}}cosI}& =-\frac{cosI{D}_{J_{\mathrm{r}}L}}{L}+\frac{D_{J_{\mathrm{r}}{L}_z}}{L}=0,\hfill \end{array}$$(D.6b)

$$\begin{array}{cc}\hfill D_{LcosI}& =-\frac{cosI{D}_{\mathit{LL}}}{L}+\frac{D_{L{L}_z}}{L}=0,\hfill \end{array}$$(D.6c)

$$\begin{array}{cc}\hfill D_{cosIcosI}& =\frac{{cos}^{2}I{D}_{\mathit{LL}}}{L^2}-\frac{2cosI{D}_{L{L}_z}}{L^2}+\frac{D_{L_z{L}_z}}{L^2},\hfill \end{array}$$(D.6d)

while the other coefficients stay unchanged. Importantly, we point out that D_{Jr cos I}​ = ​D_{L cos I}​ = ​0. This comes from the relations D_ELz​ = ​cos I D_EL and D_LLz​ = ​cos I D_LL, which are inferred from equations (B.1) and (B.5). Hence, we are left with the diffusion coefficients

$$\begin{array}{c}\hfill {\mathit{D}}_1({\mathit{J}}_{\mathrm{c}})=\left[\begin{array}{c}{D}_{J_{\mathrm{r}}}\hfill \\ D_L\hfill \\ D_{cosI}\hfill \end{array}\right];{\mathit{D}}_2({\mathit{J}}_{\mathrm{c}})=\left[\begin{array}{ccc}{D}_{J_{\mathrm{r}}{J}_{\mathrm{r}}}& D_{J_{\mathrm{r}}L}& 0\\ D_{J_{\mathrm{r}}L}& D_{\mathit{LL}}& 0\\ 0& 0& D_{cosIcosI}\end{array}\right].\end{array}$$(D.7)

Overall, they drive the FP evolution of F(J_c)​ = ​LF_rot(J), the cluster’s DF in (J_r, L, cos I).

Appendix E: N-body simulations

The simulations presented throughout the main text were performed using the direct N-body code NBODY6++GPU (Wang et al. 2015), version 4.1. We used the exact same run parameters as in appendix G of T22, and drew the initial conditions using PLUMMERPLUS.PY (Breen et al. 2017). We also refer to the aforementioned appendix for the N-body measurements, and to Table E.1 for all our binning parameters.

Each N-body realisation was composed of N​ = ​10⁵ (resp. N​ = ​10⁴) stars and integrated up to t_max​ = ​1000 HU (resp. t_max​ = ​4000 HU) with a dump every 1HU. On a 40-core CPU node with a single V100 GPU, one simulation typically required ∼​24 h (resp. ∼​12 h) of computation. Ensemble averages were performed over 50 independent runs.

In practice, the main difficulties in the N-body measurements are (i) the estimation of the instantaneous potential – necessary to compute the instantaneous radial action, J_r; (ii) the determination of the bins’ size in action space; (iii) the estimation of the relaxation rate ∂F/∂t via finite differences. This is especially true for highly tangentially anisotropic clusters (e.g., q​ = ​ − 6), where stars are closely stacked near the J_r​ = ​0 axis. For these three difficulties, we use the same approach as in T22.

Appendix F: Sphericity

In order to track the clusters’ sphericity in Fig. 2, we introduce the 3D inertia-like matrix

$$\begin{array}{c}\hfill \mathcal{I}={\displaystyle \sum _{k=1}^N}{\rho }_k^2({\mathit{r}}_k^{\mathrm{T}}{\mathit{r}}_k\mathbf{I}-{\mathit{r}}_k{\mathit{r}}_k^{\mathrm{T}})/{\displaystyle \sum _{k=1}^N}{\rho }_k^2,\end{array}$$(F.1)

with r_k the location of the k-th particle, ρ_k its local density (see Casertano & Hut 1985) and I the 3D identity matrix. In that definition, the extra ${\rho }_k^2$ factors enhance the contributions from the regions close to the centre.

The matrix ℐ is obviously symmetric. It is also semi-definite positive since for any y ∈ ℝ³, we have

$$\begin{array}{cc}\hfill {\mathbf{y}}^{\mathrm{T}}\mathcal{I}\mathbf{y}& ={\displaystyle \sum _k}{\rho }_k^2({\mathit{r}}_k^{\mathrm{T}}{\mathit{r}}_k{\mathbf{y}}^{\mathrm{T}}\mathbf{y}-{\mathbf{y}}^{\mathrm{T}}{\mathit{r}}_k{\mathit{r}}_k^{\mathrm{T}}\mathbf{y})/{\displaystyle \sum _k}{\rho }_k^2\hfill \\ \hfill & ={\displaystyle \sum _k}{\rho }_k^2(|{\mathit{r}}_k{|}^2|{\mathbf{y}}_2{|}^2-{|{\mathit{r}}_k·\mathbf{y}|}^2)/{\displaystyle \sum _k}{\rho }_k^2\ge 0,\hfill \end{array}$$(F.2)

following Cauchy–Schwarz’s inequality. As a result, ℐ has three positive eigenvalues, {λ_i}_i, which encapsulate the cluster’s sphericity. Indeed, spherically symmetric clusters have all their eigenvalues equal.

We generically define the cluster’s sphericity via h​ = ​λ_min/λ_max, which we estimated from N-body simulations. To reduce shot noise in that measurement, we averaged h over realisations as follows. First, we computed the elementary symmetric polynomials α​ = ​λ₁​ + ​λ₂​ + ​λ₃, β​ = ​λ₁λ₂​ + ​λ₁λ₃​ + ​λ₂λ₃ and γ​ = ​λ₁λ₂λ₃ for every cluster and every timestep. We then averaged the values of (α, β, γ) over all realisations. From these values, we estimated the eigenvalues λ_i as the three (positive) roots of the polynomial λ³​ − ​⟨α⟩λ²​ + ​⟨β⟩λ​ − ​⟨γ⟩.⁴

In Fig. 2, we illustrate the evolution of the sphericity, h, for various amounts of anisotropy and rotation. We typically find h​ ≃ ​0.996, that is, clusters remain reasonably spherically symmetric throughout their evolution. As as visual check of this conclusion, we represent in Fig. F.1 the late time stellar distribution of a rotating isotropic globular cluster with N​ = ​10⁴. Even close to core collapse, the cluster remains spherically symmetric.

Fig. F.1. Snapshots of the distribution of stars in a rotating (α​ = ​0.25) isotropic (q​ = ​0) Plummer cluster with N​ = ​104 stars, projected on the rotation plane. From left to right, this corresponds to t​ = ​0, 1000, 2000, 3000 HU. As time evolves, the cluster contracts and the central density increases. Notwithstanding, the cluster retains its initial spherical symmetry.

Appendix G: Rotating King model

In section 3.2, we mentioned the gravo-gyro catastrophe, introduced in Hachisu (1979) as the rotational counterpart of the gravothermal catastrophe (Lynden-Bell & Wood 1968) to explain an apparent acceleration of core collapse with rotation. To study this phenomenon, Einsel & Spurzem (1999) considered the secular evolution of a rotating King model with the initial DF

$$\begin{array}{c}\hfill F_{\mathrm{rot}}(E,L_z)\propto ({\mathrm{e}}^{-\beta E}-1){\mathrm{e}}^{-\beta {\mathrm{\Omega }}_0{L}_z},\end{array}$$(G.1)

where F_rot​ = ​F_rot(r, v)​ = ​F_rot(E, L_z) is the DF in (r, v). In this DF, $\beta \text{}=\text{}1/{\sigma }_{\text{c}}^2$ is an inverse temperature with σ_c the central velocity dispersion, while ${\omega }_0=\sqrt{9/(4\pi G{n}_{\mathrm{c}})}{\mathrm{\Omega }}_0$ is a rotation parameter with n_c the central density and Ω₀ the angular velocity in the cluster’s centre (Lagoute & Longaretti 1996). In particular, β is related to the so-called King parameter (King 1966), defined by W₀​ = ​ − β(ψ − ψ_t) where ψ_t is the potential at the cluster’s edge. The rotating King models from Einsel & Spurzem (1999) can therefore be parametrised by (W₀, ω₀). In practice, Einsel & Spurzem (1999) integrated these models up to core collapse using a FP scheme in (E, L_z) for different rotation parameters, ω₀. This is presented in their fig. 2. In particular, they demonstrated a clear correlation between the time of core collapse and the rotation parameter of the King model. The faster a cluster rotates, the shorter the core collapse time. Importantly, as emphasised in section 3.2, we recall that as one varies ω₀, the cluster’s mean potential also varies.

Appendix H: 2D Fokker-Planck equations

In practice, it is convenient to study long-term relaxation through two-dimensional projections of action space. To do so, we must integrate over one coordinate the 3D FP equation expressed either in (J_r, L, L_z) (appendix D.2) or in (J_r, L, cos I) (appendix D.3). In this appendix, we compute the 2D FP equations in (J_r, L), (J_r, cos I) and (L, cos I) that respectively drive the evolution of the DFs

$$\begin{array}{cc}\hfill F(J_{\mathrm{r}},L)& ={\displaystyle {\int }_{-L}^L}\mathrm{d}{L}_z{F}_{\mathrm{rot}}(J_{\mathrm{r}},L,L_z),\hfill \end{array}$$(H.1a)

$$\begin{array}{cc}\hfill F(J_{\mathrm{r}},cosI)& ={\displaystyle {\int }_0^{+\infty }}\mathrm{d}LF(J_{\mathrm{r}},L,cosI),\hfill \end{array}$$(H.1b)

$$\begin{array}{cc}\hfill F(L,cosI)& ={\displaystyle {\int }_0^{+\infty }}\mathrm{d}{J}_{\mathrm{r}}F(J_{\mathrm{r}},L,cosI).\hfill \end{array}$$(H.1c)

H.1. Equation in (J_r, L)

We start from equation (1), which describes diffusion in (J_r, L, L_z), and write the flux divergence as

$$\begin{array}{cc}\hfill \frac{\partial F_{\mathrm{rot}}}{\partial t}& =-(\frac{\partial {\mathcal{F}}_{J_{\mathrm{r}}}}{\partial J_{\mathrm{r}}}+\frac{\partial {\mathcal{F}}_L}{\partial L}+\frac{\partial {\mathcal{F}}_{L_z}}{\partial L_z}).\hfill \end{array}$$(H.2)

Integrating over L_z yields

$$\begin{array}{cc}\hfill {\displaystyle {\int }_{-L}^L}\mathrm{d}{L}_z\frac{\partial {\mathcal{F}}_{J_{\mathrm{r}}}}{\partial J_{\mathrm{r}}}& =\frac{\partial }{\partial J_{\mathrm{r}}}{\displaystyle {\int }_{-L}^L}\mathrm{d}{L}_z{\mathcal{F}}_{J_{\mathrm{r}}},\hfill \end{array}$$(H.3a)