11institutetext:  Institut d’Astrophysique de Paris, CNRS and Sorbonne
Université, UMR 7095, 98 bis Boulevard Arago, F-75014 Paris, France
22institutetext: Department of Physics and Astronomy, The University of North
Carolina at Chapel Hill, Chapel Hill, NC 27599, USA 33institutetext: IPhT,
DRF-INP, UMR 3680, CEA, L’Orme des Merisiers, Bât 774, 91191 Gif-sur-Yvette,
France 44institutetext: Korea Institute for Advanced Study, 85 Hoegi-ro,
Dongdaemun-gu, Seoul 02455, Republic of Korea

# Non-resonant relaxation of rotating globular clusters

Kerwann Tep E-mail<EMAIL_ADDRESS>0000-0002-0003-1873 Jean-Baptiste Fouvry 11
0000-0002-0030-371X Christophe Pichon 113344 0000-0003-0695-6735

The long-term relaxation of rotating globular clusters is investigated through
an extension of the orbit averaged Chandrasekhar non-resonant formalism. A
comparison is made with the long-term evolution of the distribution function
in action space, measured from averages of sets of $N$-body simulations up to
core collapse. The impact of rotation on in-plane relaxation is found to be
weak. In addition, we observe a clear match between theoretical predictions
and $N$-body measurements. For the class of rotating models considered, we
find no strong gravo-gyro catastrophe accelerating core collapse. Both kinetic
theory and simulations predict a reshuffling of orbital inclinations from
overpopulated regions to underpopulated ones. This trend accelerates as the
amount of rotation is increased. Yet, for orbits closer to the rotational
plane, the non-resonant prediction does not reproduce numerical measurements.
We argue that this mismatch stems from these orbits’ coherent interactions,
which are not captured by the non-resonant formalism that only addresses local
deflections.

###### Key Words.:

Diffusion – Gravitation – Galaxies: kinematics and dynamics

##  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 limited to isotropic, non-rotating, old globular
clusters with a single stellar population (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 – e.g., 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).

However, the larger scale structures in which stellar systems evolve exhibit
some degree of anisotropy inherited from the primordial density field (see,
e.g., West, 1994). These structures funnel gas infall and spin up stellar
systems either through their relative alignments (Plionis & Basilakos, 2002)
or through the infall of satellites onto the main cluster (Tormen, 1997;
Aubert et al., 2004; Pichon et al., 2011). Angular momentum may also be
generated through ejection from an initially non-rotating cold self-
gravitating systems in some cases (Benhaiem et al., 2016).

The last decade has seen the extraction of new data – e.g., HST (Bellini et
al., 2017) and Gaia DR2 (Bianchini et al., 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 which describes the evolution of globular clusters
and accounts for their rotation is needed to describe their long-term
evolution. 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, i.e. the presence of
rotation, can have a significant impact on the cluster’s long-term evolution,
e.g., 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 (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 w.r.t. 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 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
${(\bm{r},\bm{v})}$ space,
${F_{\mathrm{rot}}\\!=\\!F_{\mathrm{rot}}(\bm{r},\bm{v})}$, with $\bm{r}$ the
position and $\bm{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 in rotation. Since we are
interested in orbital distortion, it is convenient to monitor the cluster’s
evolution in action space, via
${F_{\mathrm{rot}}\\!=\\!F_{\mathrm{rot}}(\bm{J})}$, where
$\bm{J}=(J_{\mathrm{r}},L,L_{z})$ are the specific action coordinates. Here,
$J_{\mathrm{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 Fokker–Planck
(FP) equation (see, e.g., §7.4 of Binney & Tremaine, 2008)

$\displaystyle\frac{\partial F_{\mathrm{rot}}(\bm{J},t)}{\partial t}{}$
$\displaystyle\\!=\\!-\frac{\partial}{\partial\bm{J}}\\!\cdot\\!\bm{\mathcal{F}}(\bm{J})$
(1)
$\displaystyle\\!=\\!-\frac{\partial}{\partial\bm{J}}\\!\cdot\\!\bigg{[}\bm{D}_{1}(\bm{J})\,F_{\mathrm{rot}}(\bm{J})\\!-\\!\frac{1}{2}\frac{\partial}{\partial\bm{J}}\\!\cdot\\!\bigg{(}\bm{D}_{2}(\bm{J})\,F_{\mathrm{rot}}(\bm{J})\bigg{)}\bigg{]},$

with $\bm{\mathcal{F}}(\bm{J})$ the action space flux. Therein, we find the
first-order diffusion coefficient, $\bm{D}_{1}$, and the diffusion matrix,
$\bm{D}_{2}$. 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 $\bm{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

$\displaystyle\langle\Delta v_{i}\rangle\\!$
$\displaystyle=\\!-2A\\!\int\\!{\mathrm{d}}\bm{v}^{\prime}\,\frac{w_{i}}{w^{3}}F_{\mathrm{rot}}(\bm{r},\bm{v}^{\prime}),$
(2a) $\displaystyle\langle\Delta v_{i}\Delta v_{j}\rangle\\!$
$\displaystyle=\\!A\\!\int\\!{\mathrm{d}}\bm{v}^{\prime}\,\frac{w^{2}\delta_{ij}-w_{i}w_{j}}{w^{3}}F_{\mathrm{rot}}(\bm{r},\bm{v}^{\prime}),$
(2b)

where ${1\\!\leq\\!i,j\\!\leq\\!3}$ run over the three directions of the
coordinate system. We also defined ${A\\!=\\!4\pi mG^{2}\\!\ln\Lambda}$ and
the relative velocity ${\bm{w}\\!=\\!\bm{v}-\bm{v}^{\prime}}$. Here,
$\ln\Lambda$ stands for the Coulomb logarithm, set to ${\Lambda\\!=\\!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{bmatrix}\langle\Delta v_{\parallel}\rangle\\\ \langle(\Delta
v_{\parallel})^{2}\rangle\\\ \langle(\Delta
v_{\perp})^{2}\rangle\end{bmatrix}\\!=\\!A\\!\\!\int\\!\\!{\mathrm{d}}w{\mathrm{d}}\vartheta{\mathrm{d}}\phi\sin\vartheta\begin{bmatrix}-2\cos\vartheta\\\
w\sin^{2}\vartheta\\\ w\,(1+\cos^{2}\vartheta)\end{bmatrix}F_{\mathrm{rot}},$
(3)

where $\Delta v_{\parallel}$ (resp. $\Delta v_{\perp}$) is the local velocity
deflection along (resp. perpendicular to) the test star’s trajectory and we
have used polar coordinates $(w,\vartheta,\phi)$ with the $z$-axis parallel to
$\bm{v}$ (see Fig. 11). Appendix A.1 also details the arguments at which the
DF,
${F_{\mathrm{rot}}(\bm{r},\bm{v}^{\prime})\\!=\\!F_{\mathrm{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_{\mathrm{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\\!\in\\![r_{\mathrm{p}},r_{\mathrm{a}}]}$) and its angular
phase ($\varphi\\!\in\\![0,2\pi]$) within the orbital plane. This is
illustrated in Fig. 1.

Figure 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\\!\in\\![r_{\mathrm{p}},r_{\mathrm{a}}]}$) and its angular range
($\varphi\\!\in\\![0,2\pi]$), weighted by the surface density
$\Omega_{r}/(2\pi^{2}|v_{\mathrm{r}}|)$. Here, the frame is that of the right
panel of Fig. 12.

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

$D_{X}=\frac{\Omega_{r}}{\pi}\int_{r_{\mathrm{p}}}^{r_{\mathrm{a}}}\frac{{\mathrm{d}}r}{|v_{\mathrm{r}}|}\int_{0}^{2\pi}\frac{{\mathrm{d}}\varphi}{2\pi}\,\langle\Delta
X\rangle(r,\varphi),$ (4)

where $X$ runs over all the diffusion coefficients to compute, i.e. runs over
singlets and pairs from ${\\{E,L,L_{z}\\}}$. Here, $\Omega_{r}$ is the radial
frequency and $v_{\mathrm{r}}$ the radial velocity that is evaluated along the
mean field orbit. We refer to appendix C for detailed expressions. In
particular, in section C.1, we show how the angular average can be performed
explicitly for some classes of rotating models. And, in section 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
${\bm{J}\\!=\\!(J_{\mathrm{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
${\bm{J}_{\mathrm{c}}\\!=\\!(J_{\mathrm{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_{\mathrm{r}},L)}$
space, follows the same DFs 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

$F_{\mathrm{rot}}(\bm{J})=F_{\mathrm{tot}}(J_{\mathrm{r}},L)\,\big{(}1+\alpha\,\mathrm{sgn}[L_{z}/L]\big{)},$
(5)

where $F_{\mathrm{tot}}(J_{\mathrm{r}},L)$ is the non-rotating DF and $\alpha$
a dimensionless parameter between 0 and 1. Physically, this parameter
corresponds to converting a fraction $\alpha$ of retrograde orbits
${(L_{z}\\!<\\!0)}$ into prograde orbits ${(L_{z}\\!>\\!0)}$. Importantly,
since ${L_{z}\\!\mapsto\\!\mathrm{sgn}[L_{z}/L]}$ is an odd function, we
stress that $F_{\mathrm{tot}}$ and $F_{\mathrm{rot}}$ generate the exact same
potential.

When considering orbital inclinations via
${\bm{J}_{\mathrm{c}}\\!=\\!(J_{\mathrm{r}},L,\cos I)}$ (see appendix D.3),
the LBD yields the DF in $\bm{J}_{\mathrm{c}}$ space 111Using $\cos I$ as an
effective coordinate to describe relaxation might be problematic for some
clusters. We refer to Appendix J for more details.

$F(\bm{J}_{\mathrm{c}})=L\,F_{\mathrm{tot}}(J_{\mathrm{r}},L)\,\big{(}1\\!+\\!\alpha\,\mathrm{sgn}[\cos
I]\big{)}.$ (6)

Integrating this equation over $\cos I$, one recovers the reduced DF in
${(J_{\mathrm{r}},L)}$,
${F(J_{\mathrm{r}},L)\\!=\\!2LF_{\mathrm{tot}}(J_{\mathrm{r}},L)}$, used to
describe non-rotating clusters (see, e.g., Hamilton et al., 2018).

###  3.1 Sphericality 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 restrict ourselves to (linearly)
stable rotating clusters, i.e. clusters that would remain spherically
symmetric.

To probe the conservation of spherical symmetry, we investigate the
sphericality, $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.

Figure 2: Sphericality, $h$, of a sample of rotating anisotropic clusters, as
defined in appendix F, from $N$-body simulations for ${N\\!=\\!10^{5}}$. 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
${\alpha\\!=\\!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,\alpha)\\!=\\!(1,0.5)}$, and in a smaller fashion
${(q,\alpha)\\!=\\!(1,0.25)}$ – while the others are linearly stable and
remain spherically symmetric.

Except for the radially-anisotropic and rapidly-rotating cluster
${(q,\alpha)\\!=\\!(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,\alpha)$ 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

$t_{\mathrm{rh}}=\frac{0.138\,N\,r_{\mathrm{h}}^{3/2}}{(GM)^{1/2}\ln(0.11\,N)},$
(7)

with $r_{\mathrm{h}}$ the half-mass radius. For a cluster with
${N\\!=\\!10^{5}}$ stars, this yields
${t_{\mathrm{rh}}\\!=\\!994\,{\mathrm{HU}}}$, where ${\mathrm{HU}}$ stands for
the Hénon units (Heggie & Mathieu, 1986) in which
${G\\!=\\!M\\!=\\!R_{\mathrm{v}}\\!=\\!1}$, with $R_{\mathrm{v}}$ the virial
radius222For a Plummer cluster, this is related to the Plummer scale length
$b$ by the relation ${R_{\mathrm{v}}\\!=\\!16b/(3\pi)}$ (see, e.g., Dejonghe,
1987).. In practice, with ${N\\!=\\!10^{5}}$ running a numerical simulation up
to ${t\\!\sim\\!t_{\mathrm{rh}}}$ takes about one day of computation, see
appendix E. In Breen et al. (2017), core collapse occurs at
${t\\!\sim\\!\\!17t_{\mathrm{rh}}}$ for an isotropic cluster (see table 1
therein for the dependence w.r.t. velocity anisotropy). Hence, integrating
clusters with ${N\\!=\\!10^{5}}$ up to core collapse is not reasonably
feasible. In this section we therefore scale down the size of the clusters to
${N\\!=\\!10^{4}}$ stars. In that case,
${t_{\mathrm{rh}}\\!=\\!132\,{\mathrm{HU}}}$, and core collapse is numerically
reached in about ten hours. Since we will focus on measuring the core radius –
a very integrated quantity – the quality of the numerical measurements will
not be too much degraded by this use of a smaller value for $N$. In practice,
measurements are averaged over 50 runs.

Figure 3 represents the time evolution of the averaged core radius (see
equations II.2–II.4 of Casertano & Hut, 1985) for the usual sets of
anisotropies and rotations (${q\\!=\\!1,0,-6}$ and
${\alpha\\!=\\!0,0.1,0.25,0.5}$).

Figure 3: Core radius as a function of time, as measured in $N$-body
simulations, with $t_{\mathrm{rh}}$ defined in equation (7). Each case has
been averaged over 50 realisations with ${N\\!=\\!10^{4}}$. Increasing the
rotation strength $\alpha$ 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).

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_{0}\\!=\\!6.0}$, a non-rotating King cluster collapses at
${t\sim 12t_{\mathrm{rh}}}$, while a rotating one, with
${\omega_{0}/\Omega_{0}=0.4}$, can collapse as early as ${t\sim
9t_{\mathrm{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 $\omega_{0}$ (King) and $\alpha$
(Plummer) do not parametrise rotation in the same way. As Einsel & Spurzem
(1999) varies $\omega_{0}$, the mean density profile gets modified. Indeed,
the rotating DF of the King model (equation 73) 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 J). This is left for future investigations.

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

Using once again $N$-body simulations with ${N\\!=\\!10^{4}}$ stars, let us
finally investigate the typical features of the in-plane diffusion, i.e.
relaxation in ${(J_{\mathrm{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$, i.e. the in-
plane relaxation of the distribution of angular momenta.

Figure 4: Evolution of the angular momentum DF, $F(L)$, as measured in
$N$-body simulations with ${N\\!=\\!10^{4}}$, with $t_{\mathrm{rh}}$ 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.

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.

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$.

Figure 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.

As could have been expected, the clusters tend towards smoother distribution
of inclinations, i.e. relaxation tends to erase discontinuities. The sharp
discontinuity at ${\cos I\\!=\\!0}$ has already been washed away after
$1\,t_{\mathrm{rh}}$. Additionally, the larger the initial rotation, the
faster the disappearance of the discontinuity. 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.

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 $\bm{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” i.e. relaxation
occurring in ${(J_{\mathrm{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^{5}}$ stars.

Figure 6: Initial evolution of the core radius in clusters with
${N\\!=\\!10^{5}}$ 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.

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, e.g., for ${(q,\alpha)\\!=\\!(1,0.5)}$ (see Fig. 2), we find
numerically (not reported in Fig. 6) that the clusters flatten and that
relaxation is accelerated.

To probe more precisely relaxation, let us now measure the relaxation rate,
${\partial F/\partial t}$, in ${(J_{\mathrm{r}},L)}$. This is presented in
Fig. 7 with details spelled out in Appendix E.

Figure 7: Relaxation rate, $\partial F/\partial t$, in $(J_{\mathrm{r}},L)$
for various velocity anisotropies $q$ (left to right) and rotation parameters
$\alpha$ (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 $\alpha$.

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, $\alpha$, has no significant impact on
the geometry of the in-plane diffusion. Increasing $\alpha$ only (very) weakly
accelerates the in-plane relaxation. This is in concordance with the slightly
shorter core collapse time observed in Fig. 3.

###  4.2 Non-resonant prediction

Let us now compute the NR prediction for the relaxation rate in
${(J_{\mathrm{r}},L)}$. This prediction is numerically very intensive and
requires 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) allows us to analytically perform the angular part of the orbit
average (equation 4). This helps reducing the computation time. We refer to
Appendix C.1 for details.

The relaxation rate, ${\partial F/\partial t}$ in ${(J_{\mathrm{r}},L)}$,
presented in Fig. 8 for various anisotropies and rotations,

Figure 8: Relaxation rate, $\partial F/\partial t$, in ${(J_{\mathrm{r}},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 which depends on $q$, and
(very) weakly on $\alpha$.

suggests once again an orbital reshuffling toward isotropisation. In addition,
increasing the rotation parameter $\alpha$ only very weakly increases the
relaxation rate in $(J_{\mathrm{r}},L)$, while barely impacting the structures
observed in the non-rotating case.

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 which 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_{\mathrm{r}},\cos I)}$, and we refer
to appendix I 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.

Figure 9: Relaxation rate, $\partial F/\partial t$, in ${(J_{\mathrm{r}},\cos
I)}$ for various anisotropies $q$ (left to right) and rotation $\alpha$ (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 expected, in the absence of rotation (top row), the relaxation is
independent of $\cos I$, and the distribution of orientations remains uniform.

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 $\alpha$, 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_{\mathrm{r}}$ is observed, whatever $\cos I$.
Conversely, the number of orbits systematically increases for large
$J_{\mathrm{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_{\mathrm{r}}$. Overall, all clusters, independently of their anisotropies,
evolve towards smoother distribution of orbital inclinations.

As a complement of Fig. 9, appendix I also presents the relaxation rates in
${(L,\cos I)}$. It reaches similar results as for diffusion in
${(J_{\mathrm{r}},\cos I)}$. Namely, the in-plane distributions tend to
isotropise and orbital inclinations diffuse so as to reach smoother
distributions.

###  5.2 Non-resonant prediction

Let us now perform the same investigation in ${(J_{\mathrm{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.

Figure 10: Relaxation rate, $\partial F/\partial t$, in ${(J_{\mathrm{r}},\cos
I)}$ as in Fig. 10, 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.

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. 10) 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, i.e. as one considers ${\cos
I\\!\to\\!\pm 1}$. This is in sharp contrast with the “straighter” structures
predicted by the NR theory (Fig. 9), where the relaxation rate does not vary
so much as a function of orbital inclination. This seems to be a robust
conclusion, as this mismatch remains as one increases the sampling of the
various integrations performed in the NR prediction. A similar conclusion is
reached in appendix I 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 ${\alpha\\!\neq\\!0}$, recall that ${\cos
I\\!=\\!\pm 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,
${\langle\bm{L}\rangle\\!\propto\\!\alpha\,\mathbf{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. This is
likely why we obtain a prediction which only weakly depends on orbital
inclination (Fig. 10).

Following the work of Meiron & Kocsis (2019) (see also references therein), we
anticipate for VRR to play a key role in these coherent interactions. Indeed,
VRR, which is driven by persistent torques between orbital planes, might be
needed to explain (i) the faster azimuthal 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 predict 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, i.e. 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. However, we noted that for for orbits aligned with the cluster’s
rotation plane, 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, e.g. as demonstrated 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, i.e. 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. Indeed, as argued in section 5, rotation reduces the
respective phase mixing of orbital planes, hence increasing the strength of
collective effects. 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).

Eventually, understanding statistically (beyond the mean response, see
Touchette, 2009; Feliachi & Fouvry, 2023) the impact of rotation on secular
evolution should prove critical to model the morphological evolution of
populations of galaxies across cosmic times (see, e.g., Pichon & Aubert, 2006;
Conselice, 2014). Indeed, the distribution of angular momentum is, beyond
mass, the most relevant dynamical parameter for such populations: it separates
elliptical galaxies from spiral ones. The non-linear long-term impact of
rotation on galactic orbital structures should be the focus of future
research.

## Data Distribution

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.

## Acknowledgements

This work is partially supported by the grant _SEGAL_ 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 ‘ _CosmicWeb23_ ’ during which
this project was advanced. This work is partially supported by the National
Science Foundation under Grant No. NSF PHY-1748958. 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.

## References

  * Aarseth et al. (1974) Aarseth, S. J., Hénon, M., & Wielen, R. 1974, A&A, 37, 183
  * Akiyama & Sugimoto (1989) Akiyama, K. & Sugimoto, D. 1989, PASJ, 41, 991
  * Aubert et al. (2004) Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376
  * Bar-Or & Alexander (2016) Bar-Or, B. & Alexander, T. 2016, ApJ, 820, 129
  * Bellini et al. (2017) Bellini, A. et al. 2017, ApJ, 842, 6
  * Benhaiem et al. (2016) Benhaiem, D., Joyce, M., Sylos Labini, F., & Worrakitpoonpon, T. 2016, A&A, 585, A139
  * Bezanson et al. (2017) Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. 2017, SIAM Review, 59, 65
  * Bianchini et al. (2013) Bianchini, P., Varri, A. L., Bertin, G., & Zocchi, A. 2013, ApJ, 772, 67
  * Bianchini et al. (2018) Bianchini, P. et al. 2018, MNRAS, 481, 2125
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton Univ. Press)
  * Breen et al. (2017) Breen, P. G., Varri, A. L., & Heggie, D. C. 2017, MNRAS, 471, 2778
  * Casertano & Hut (1985) Casertano, S. & Hut, P. 1985, ApJ, 298, 80
  * Chandrasekhar (1943) Chandrasekhar, S. 1943, ApJ, 97, 255
  * Cohn (1979) Cohn, H. 1979, ApJ, 234, 1036
  * Conselice (2014) Conselice, C. J. 2014, ARA&A, 52, 291
  * Dejonghe (1986) Dejonghe, H. 1986, Phys. Rep., 133, 217
  * Dejonghe (1987) Dejonghe, H. 1987, MNRAS, 224, 13
  * Dejonghe & de Zeeuw (1988) Dejonghe, H. & de Zeeuw, T. 1988, ApJ, 333, 90
  * Dekel et al. (2023) Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, 523, 3201
  * Eilon et al. (2009) Eilon, E., Kupi, G., & Alexander, T. 2009, ApJ, 698, 641
  * Einsel & Spurzem (1999) Einsel, C. & Spurzem, R. 1999, MNRAS, 302, 81
  * Ernst et al. (2007) Ernst, A., Glaschke, P., Fiestas, J., Just, A., & Spurzem, R. 2007, MNRAS, 377, 465
  * Fabricius et al. (2014) Fabricius, M. H., Noyola, E., Rukdee, S., et al. 2014, ApJL, 787, L26
  * Feliachi & Fouvry (2023) Feliachi, O. & Fouvry, J.-B. 2023, arXiv, 2308.08308
  * Ferraro et al. (2018) Ferraro, F. R., Mucciarelli, A., Lanzoni, B., et al. 2018, ApJ, 860, 50
  * Fiestas & Spurzem (2010) Fiestas, J. & Spurzem, R. 2010, MNRAS, 405, 194
  * Foote et al. (2020) Foote, H. R., Generozov, A., & Madigan, A.-M. 2020, ApJ, 890, 175
  * Fouvry et al. (2019) Fouvry, J.-B., Bar-Or, B., & Chavanis, P.-H. 2019, ApJ, 883, 161
  * Fouvry et al. (2021) Fouvry, J.-B., Hamilton, C., Rozier, S., & Pichon, C. 2021, MNRAS, 508, 2210
  * Geyer et al. (1983) Geyer, E. H., Hopp, U., & Nelles, B. 1983, A&A, 125, 359
  * Giersz & Heggie (1994) Giersz, M. & Heggie, D. C. 1994, MNRAS, 268, 257
  * Ginat et al. (2023) Ginat, Y. B., Panamarev, T., Kocsis, B., & Perets, H. B. 2023, MNRAS, 525, 4202
  * Greene & Ho (2004) Greene, J. E. & Ho, L. C. 2004, ApJ, 610, 722
  * Greene et al. (2020) Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257
  * Gruzinov et al. (2020) Gruzinov, A., Levin, Y., & Zhu, J. 2020, ApJ, 905, 11
  * Hachisu (1979) Hachisu, I. 1979, PASJ, 31, 523
  * Hachisu (1982) Hachisu, I. 1982, PASJ, 34, 313
  * Hamilton et al. (2018) Hamilton, C., Fouvry, J.-B., Binney, J., & Pichon, C. 2018, MNRAS, 481, 2041
  * Heggie & Hut (2003) Heggie, D. & Hut, P. 2003, The Gravitational Million-Body Problem (Cambridge University Press)
  * Heggie & Mathieu (1986) Heggie, D. C. & Mathieu, R. D. 1986, in The Use of Supercomputers in Stellar Dynamics, ed. P. Hut & S. L. W. McMillan, Vol. 267 (Springer Berlin), 233
  * Heyvaerts (2010) Heyvaerts, J. 2010, MNRAS, 407, 355
  * Inagaki & Hachisu (1978) Inagaki, S. & Hachisu, I. 1978, PASJ, 30, 39
  * Jindal et al. (2019) Jindal, A., Webb, J. J., & Bovy, J. 2019, MNRAS, 487, 3693
  * Kamann et al. (2018) Kamann, S., Husser, T.-O., Dreizler, S., et al. 2018, MNRAS, 473, 5591
  * Kamlah et al. (2022) Kamlah, A. W. H., Spurzem, R., Berczik, P., et al. 2022, MNRAS, 516, 3266
  * Kim et al. (2004) Kim, E., Lee, H. M., & Spurzem, R. 2004, MNRAS, 351, 220
  * King (1966) King, I. R. 1966, AJ, 71, 64
  * Kocsis & Tremaine (2011) Kocsis, B. & Tremaine, S. 2011, MNRAS, 412, 187
  * Kocsis & Tremaine (2015) Kocsis, B. & Tremaine, S. 2015, MNRAS, 448, 3265
  * Kontizas et al. (1989) Kontizas, E., Kontizas, M., Sedmak, G., & Smareglia, R. 1989, AJ, 98, 590
  * Kormendy & Ho (2013) Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511
  * Lagoute & Longaretti (1996) Lagoute, C. & Longaretti, P. Y. 1996, A&A, 308, 441
  * Lanzoni et al. (2018) Lanzoni, B., Ferraro, F. R., Mucciarelli, A., et al. 2018, ApJ, 865, 11
  * Livernois et al. (2022) Livernois, A. R., Vesperini, E., Varri, A. L., Hong, J., & Tiongco, M. 2022, MNRAS, 512, 2584
  * Lynden-Bell (1960) Lynden-Bell, D. 1960, MNRAS, 120, 204
  * Lynden-Bell & Wood (1968) Lynden-Bell, D. & Wood, R. 1968, MNRAS, 138, 495
  * Magnan et al. (2022) Magnan, N., Fouvry, J.-B., Pichon, C., & Chavanis, P.-H. 2022, MNRAS, 514, 3452
  * McLaughlin & van der Marel (2005) McLaughlin, D. E. & van der Marel, R. P. 2005, ApJs, 161, 304
  * Meiron & Kocsis (2019) Meiron, Y. & Kocsis, B. 2019, ApJ, 878, 138
  * Miocchi et al. (2013) Miocchi, P., Lanzoni, B., Ferraro, F. R., et al. 2013, ApJ, 774, 151
  * Pichon & Aubert (2006) Pichon, C. & Aubert, D. 2006, MNRAS, 368, 1657
  * Pichon et al. (2011) Pichon, C., Pogosyan, D., Kimm, T., et al. 2011, MNRAS, 418, 2493
  * Plionis & Basilakos (2002) Plionis, M. & Basilakos, S. 2002, MNRAS, 329, L47
  * Rauch & Tremaine (1996) Rauch, K. P. & Tremaine, S. 1996, New. Astron., 1, 149
  * Risken (1996) Risken, H. 1996, The Fokker-Planck Equation, 2nd edn. (Berlin: Springer), gst:
  * Rosenbluth et al. (1957) Rosenbluth, M. N., MacDonald, W. M., & Judd, D. L. 1957, Phys. Rev., 107, 1
  * Rozier et al. (2019) Rozier, S., Fouvry, J.-B., Breen, P. G., et al. 2019, MNRAS, 487, 711
  * Sollima et al. (2019) Sollima, A., Baumgardt, H., & Hilker, M. 2019, MNRAS, 485, 1460
  * Spitzer (1975) Spitzer, L. 1975, in Dynamics of the Solar Systems, Vol. 69, 3
  * Stetson et al. (2019) Stetson, P. B., Pancino, E., Zocchi, A., Sanna, N., & Monelli, M. 2019, MNRAS, 485, 3042
  * Szölgyén & Kocsis (2018) Szölgyén, Á. & Kocsis, B. 2018, PRL, 121, 101101
  * Szölgyén et al. (2021) Szölgyén, Á., Máthé, G., & Kocsis, B. 2021, ApJ, 919, 140
  * Szölgyén et al. (2019) Szölgyén, Á., Meiron, Y., & Kocsis, B. 2019, ApJ, 887, 123
  * Tep et al. (2022) Tep, K., Fouvry, J.-B., & Pichon, C. 2022, MNRAS, 514, 875
  * Tiongco et al. (2021) Tiongco, M., Collier, A., & Varri, A. L. 2021, MNRAS, 506, 4488
  * Tiongco et al. (2016) Tiongco, M., Vesperini, E., & Varri, A. L. 2016, MNRAS, 461, 402
  * Tiongco et al. (2020) Tiongco, M., Vesperini, E., & Varri, A. L. 2020, in Star Clusters: From the Milky Way to the Early Universe, Vol. 351, 524–527
  * Tiongco et al. (2022) Tiongco, M., Vesperini, E., & Varri, A. L. 2022, MNRAS, 512, 1584
  * Tormen (1997) Tormen, G. 1997, MNRAS, 290, 411
  * Touchette (2009) Touchette, H. 2009, Phys. Rep, 478, 1
  * Trager et al. (1995) Trager, S. C., King, I. R., & Djorgovski, S. 1995, AJ, 109, 218
  * Varri & Bertin (2012) Varri, A. L. & Bertin, G. 2012, A&A, 540, A94
  * Vasiliev (2015) Vasiliev, E. 2015, MNRAS, 446, 3150
  * Wang et al. (2015) Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070
  * Watkins et al. (2015) Watkins, L. L., van der Marel, R. P., Bellini, A., & Anderson, J. 2015, ApJ, 803, 29
  * West (1994) West, M. J. 1994, MNRAS, 268, 79
  * White & Shawl (1987) White, R. E. & Shawl, S. J. 1987, ApJ, 317, 246
  * Wilson (1975) Wilson, C. P. 1975, AJ, 80, 175

##  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. 11.

Figure 11: Tailored frame used to compute the parallel and perpendicular local
velocity deflections in equations (8). By construction, the test star’s
angular momentum $\bm{L}$ is along the axis $\mathbf{e}_{2}$.

The local velocity deflections are then given by

$\displaystyle\langle\Delta v_{\parallel}\rangle$ $\displaystyle=\langle\Delta
v_{1}\rangle,$ (8a) $\displaystyle\langle(\Delta v_{\parallel})^{2}\rangle$
$\displaystyle=\langle(\Delta v_{1})^{2}\rangle,$ (8b)
$\displaystyle\langle(\Delta v_{\perp})^{2}\rangle$
$\displaystyle=\langle(\Delta v_{2})^{2}\rangle+\langle(\Delta
v_{3})^{2}\rangle,$ (8c)

where the coordinate 1 is the $z$ coordinate, parallel to $\bm{v}$. The
coordinates 2 and 3 are defined such that the projection of $\bm{r}$ on the
$(Oxy)$ plane is along $\mathbf{e}_{3}$ (Fig. 11). We now introduce
${\bm{w}\\!=\\!\bm{v}\\!-\\!\bm{v}^{\prime}}$ and change the integration
variables from $\bm{v}^{\prime}$ to $\bm{w}$ in equation (2). We get

$\displaystyle\langle\Delta v_{i}\rangle$
$\displaystyle=-2A\int{\mathrm{d}}\bm{w}\,\frac{w_{i}}{w^{3}}F_{\mathrm{rot}},$
(9a) $\displaystyle\langle\Delta v_{i}\Delta v_{j}\rangle$
$\displaystyle=A\int{\mathrm{d}}\bm{w}\,\frac{w^{2}\delta_{ij}-w_{i}w_{j}}{w^{3}}F_{\mathrm{rot}},$
(9b)

with the shortened notation
${F_{\mathrm{rot}}\\!=\\!F_{\mathrm{rot}}(r,\bm{v}^{\prime}\\!=\\!\bm{v}\\!-\\!\bm{w})}$.
Using these equations allows us to compute the needed coefficients from
equation (8) through

$\begin{bmatrix}\langle\Delta v_{1}\rangle\\\ \langle(\Delta
v_{1})^{2}\rangle\\\ \langle(\Delta v_{2})^{2}\rangle\\\ \langle(\Delta
v_{3})^{2}\rangle\end{bmatrix}=A\\!\\!\int\\!\\!\frac{{\mathrm{d}}\bm{w}}{w^{3}}\begin{bmatrix}-2w_{1}\\\\[2.15277pt]
w^{2}-w_{1}^{2}\\\\[2.15277pt] w^{2}-w_{2}^{2}\\\\[2.15277pt]
w^{2}-w_{3}^{2}\end{bmatrix}F_{\mathrm{rot}}.$ (10)

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

$\displaystyle
w_{1}=w\cos\vartheta\,;\,w_{2}=w\sin\vartheta\cos\phi\,;\,w_{3}=w\sin\vartheta\sin\phi.$
(11)

Injecting equations (11) into equations (10), we obtain

$\begin{bmatrix}\langle\Delta v_{1}\rangle\\\ \langle(\Delta
v_{1})^{2}\rangle\\\ \langle(\Delta v_{2})^{2}\rangle\\\ \langle(\Delta
v_{3})^{2}\rangle\end{bmatrix}\\!=\\!A\\!\\!\int\\!\\!{\mathrm{d}}w{\mathrm{d}}\vartheta{\mathrm{d}}\phi\sin\vartheta\begin{bmatrix}-2\cos\vartheta\\\
w\sin^{2}\vartheta\\\ w(1-\sin^{2}\vartheta\cos^{2}\phi)\\\
w(1-\sin^{2}\vartheta\sin^{2}\phi)\end{bmatrix}F_{\mathrm{rot}}.$ (12)

Here, the background distribution, $F_{\mathrm{rot}}$, is to be evaluated in
${F_{\mathrm{rot}}(\bm{r},\bm{v}^{\prime})\\!=\\!F_{\mathrm{rot}}(E^{\prime},L^{\prime},L_{z}^{\prime})}$.
We compute these arguments in appendix A.2. Finally injecting equation (12)
into equation (8), 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^{\prime}$, $L^{\prime}$ and
$L_{z}^{\prime}$ for the background distribution function,
$F_{\mathrm{rot}}(E^{\prime},L^{\prime},L_{z}^{\prime})$, in equation (2).

We start with $E^{\prime}$. It reads

$\displaystyle E^{\prime}$ $\displaystyle\\!=\\!\psi(r)+\frac{\bm{v}^{\prime
2}}{2}\\!=\\!\psi(r)+\frac{(\bm{v}-\bm{w})^{2}}{2}\\!=\\!E+\frac{w^{2}}{2}-wv\cos\vartheta,$
(13)

since ${\bm{v}\\!=\\!v\,\mathbf{e}_{1}}$ (see Fig. 11). This form allows us to
get an upper bound on the $w$-integral of equation (12), above which
${E^{\prime}\\!>\\!0}$. As soon as ${E^{\prime}\\!>\\!0}$, the background star
is unbound and its DF vanishes.

The computation of $L^{\prime}$ and $L_{z}^{\prime}$ is best achieved by using
a different coordinate system, as illustrated in Figs. 11 and 12.

Figure 12: Illustration of the frames used to compute the argument of the
background DF, $F_{\mathrm{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.

More precisely, the frame $(1,2,3)$ is related to the $(X,Y,Z)$ one through
the relations

$\left[\begin{array}[]{l}\mathbf{e}_{1}\\\ \mathbf{e}_{2}\\\
\mathbf{e}_{3}\end{array}\right]=\left[\begin{array}[]{ccc}\frac{v_{\mathrm{r}}}{v}\cos\varphi-\frac{v_{\mathrm{t}}}{v}\sin\varphi&\frac{v_{\mathrm{r}}}{v}\sin\varphi+\frac{v_{\mathrm{t}}}{v}\cos\varphi&0\\\
0&0&1\\\
\frac{v_{\mathrm{r}}}{v}\sin\varphi+\frac{v_{\mathrm{t}}}{v}\cos\varphi&\frac{v_{\mathrm{t}}}{v}\sin\varphi-\frac{v_{\mathrm{r}}}{v}\cos\varphi&0\end{array}\right]\left[\begin{array}[]{l}\displaystyle{\mathbf{e}_{X}}\\\
\displaystyle{\mathbf{e}_{Y}}\\\
\displaystyle{\mathbf{e}_{Z}}\end{array}\right],$ (14)

with $v_{\mathrm{r}}$ (resp. $v_{\mathrm{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

$\displaystyle\left[\begin{array}[]{l}\mathbf{e}_{X}\\\ \mathbf{e}_{Y}\\\
\mathbf{e}_{Z}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{ccc}\cos I&0&\sin I\\\ 0&1&0\\\ -\sin
I&0&\cos
I\end{array}\right]\left[\begin{array}[]{l}\displaystyle{\mathbf{e}_{x}}\\\
\displaystyle{\mathbf{e}_{y}}\\\
\displaystyle{\mathbf{e}_{z}}\end{array}\right].$ (24)

We can now compute $L$, the norm of $\bm{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 (10). We write

$\displaystyle L^{\prime 2}$
$\displaystyle\\!=\\!(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},$
(25)

where $(r_{1},r_{2},r_{3})$ are the coordinates of $\bm{r}$ in the frame
$(1,2,3)$. Using equations (14) and (24), we obtain

$r_{1}=\frac{v_{\mathrm{r}}}{v}r\quad;\quad r_{2}=0\quad;\quad
r_{3}=\frac{v_{\mathrm{t}}}{v}r.$ (26)

Therefore, one gets

$\displaystyle L^{\prime}$ $\displaystyle=r\sqrt{v_{2}^{\prime
2}+\bigg{(}\frac{v_{\mathrm{r}}}{v}v_{3}^{\prime}-\frac{v_{\mathrm{t}}}{v}v_{1}^{\prime}\bigg{)}^{2}}$
(27)
$\displaystyle=r\sqrt{(w\sin\vartheta\cos\phi)^{2}+\bigg{(}v_{\mathrm{t}}+\frac{v_{\mathrm{r}}}{v}w\sin\vartheta\sin\phi-\frac{v_{\mathrm{t}}}{v}w\cos\vartheta\bigg{)}^{2}},$

given that $\bm{v}^{\prime}\\!=\\!\bm{v}-\bm{w}$.

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

$L_{z}^{\prime}=xv_{y}^{\prime}-yv_{x}^{\prime}.$ (28)

From equation (24), the positions read

$\displaystyle x$ $\displaystyle=X\cos I=r\cos\varphi\cos I,$ (29a)
$\displaystyle y$ $\displaystyle=Y=r\sin\varphi,$ (29b) $\displaystyle z$
$\displaystyle=X\sin I=r\cos\varphi\sin I,$ (29c)

using ${X\\!=\\!r\cos\varphi}$ and ${Y\\!=\\!r\sin\varphi}$. We then obtain
the velocities

$\displaystyle v_{x}$ $\displaystyle=v_{\mathrm{r}}\cos\varphi\cos
I-v_{\mathrm{t}}\sin\varphi\cos I,$ (30a) $\displaystyle v_{y}$
$\displaystyle=v_{\mathrm{r}}\sin\varphi+v_{\mathrm{t}}\cos\varphi,$ (30b)
$\displaystyle v_{z}$ $\displaystyle=v_{\mathrm{r}}\cos\varphi\sin
I-v_{\mathrm{t}}\sin\varphi\sin I,$ (30c)

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

$\displaystyle\left[\begin{array}[]{l}v^{\prime}_{x}\\\ v^{\prime}_{y}\\\
v^{\prime}_{z}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{ccc}\cos I&0&-\sin I\\\ 0&1&0\\\ \sin
I&0&\cos
I\end{array}\right]\left[\begin{array}[]{l}\displaystyle{v^{\prime}_{X}}\\\
\displaystyle{v^{\prime}_{Y}}\\\
\displaystyle{v^{\prime}_{Z}}\end{array}\right].$ (40)

Second, using the relation in equations (14) (Fig. 11 and right panel of Fig.
12), we have

$\left[\begin{array}[]{l}v^{\prime}_{X}\\\ v^{\prime}_{Y}\\\
v^{\prime}_{Z}\end{array}\right]=\left[\begin{array}[]{ccc}\frac{v_{\mathrm{r}}}{v}\cos\varphi-\frac{v_{\mathrm{t}}}{v}\sin\varphi&0&\frac{v_{\mathrm{r}}}{v}\sin\varphi+\frac{v_{\mathrm{t}}}{v}\cos\varphi\\\
\frac{v_{\mathrm{r}}}{v}\sin\varphi+\frac{v_{\mathrm{t}}}{v}\cos\varphi&0&\frac{v_{\mathrm{t}}}{v}\sin\varphi-\frac{v_{\mathrm{r}}}{v}\cos\varphi\\\
0&1&0\end{array}\right]\left[\begin{array}[]{l}\displaystyle{v^{\prime}_{1}}\\\
\displaystyle{v^{\prime}_{2}}\\\
\displaystyle{v^{\prime}_{3}}\end{array}\right].$ (41)

Overall, we can rewrite equation (28) as

$\displaystyle L_{z}^{\prime}$
$\displaystyle=r\,\bigg{(}\frac{v^{\prime}_{1}v_{\mathrm{t}}}{v}-\frac{v^{\prime}_{3}v_{\mathrm{r}}}{v}\bigg{)}\cos
I+rv^{\prime}_{2}\sin\varphi\sin I$
$\displaystyle=r\,\bigg{(}v_{\mathrm{t}}-\frac{v_{\mathrm{t}}}{v}w\cos\vartheta+\frac{v_{\mathrm{r}}}{v}w\sin\vartheta\sin\phi\bigg{)}\cos
I$ $\displaystyle\quad-rw\sin\vartheta\cos\phi\sin\varphi\sin I.$ (42)

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

$\displaystyle|L_{z}^{\prime}|\leq
L^{\prime}\sqrt{\cos^{2}I+\sin^{2}\varphi\sin^{2}I},$ (43)

which ensures that ${|L_{z}^{\prime}|\\!\leq\\!L^{\prime}}$.

Equations (13), (27) and (42) are the main results of this section. They
provide us with explicit expressions for the arguments at which to evaluate
${F_{\mathrm{rot}}\\!=\\!F_{\mathrm{rot}}(E^{\prime},L^{\prime},L_{z}^{\prime})}$
in equation (3).

##  Appendix B 3D diffusion coefficients

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

$\displaystyle\langle\Delta E\rangle$
$\displaystyle=\tfrac{1}{2}\langle(\Delta
v_{\parallel})^{2}\rangle+\tfrac{1}{2}\langle(\Delta
v_{\perp})^{2}\rangle+v\langle\Delta v_{\parallel}\rangle,$ (44a)
$\displaystyle\langle(\Delta E)^{2}\rangle$ $\displaystyle=v^{2}\langle(\Delta
v_{\parallel})^{2}\rangle,$ (44b) $\displaystyle\langle\Delta L\rangle$
$\displaystyle=r\frac{v_{\mathrm{t}}}{v}\langle\Delta
v_{\parallel}\rangle+\frac{r^{2}}{4L}\langle(\Delta v_{\perp})^{2}\rangle,$
(44c) $\displaystyle\langle(\Delta L^{2})\rangle$
$\displaystyle=r^{2}\frac{v_{\mathrm{t}}^{2}}{v^{2}}\langle(\Delta
v_{\parallel})^{2}\rangle+\frac{r^{2}}{2}\frac{v_{\mathrm{r}}^{2}}{v^{2}}\langle(\Delta
v_{\perp})^{2}\rangle,$ (44d) $\displaystyle\langle\Delta E\Delta L\rangle$
$\displaystyle=L\langle(\Delta v_{\parallel})^{2}\rangle.$ (44e)

In order to compute the 3D FP equation (1), we also need the local diffusion
coefficients involving ${L_{z}\\!=\\!\bm{L}\\!\cdot\\!\mathbf{e}_{z}}$ .
Performing a first-order variation of ${\bm{L}\\!=\\!\bm{r}\times\bm{v}}$, we
can write

$\displaystyle\Delta\bm{L}$ $\displaystyle=\bigg{(}\frac{\Delta
v_{\parallel}}{v}-\frac{v_{\mathrm{r}}}{v_{\mathrm{t}}}\frac{\Delta
v_{3}}{v}\bigg{)}\,\bm{L}+L\frac{\Delta
v_{2}}{v}\bigg{(}\frac{v_{\mathrm{r}}}{v_{\mathrm{t}}}\mathbf{e}_{3}-\hat{\bm{v}}\bigg{)},$
(45a) $\displaystyle\mathbf{e}_{3}$
$\displaystyle=\frac{\bm{v}\times\bm{L}}{|\bm{v}\times\bm{L}|}=\frac{v\hat{\bm{r}}-v_{\mathrm{r}}\hat{\bm{v}}}{v_{\mathrm{t}}},$
(45b)

where $(\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{e}_{3})$ is the frame
illustrated in Fig. 11. We can then write

$\mathbf{e}_{3}\cdot\mathbf{e}_{z}=\sin
I\,\bigg{(}\frac{v_{\mathrm{r}}}{v}\sin\varphi+\frac{v_{\mathrm{t}}}{v}\cos\varphi\bigg{)}.$
(46)

As a result, it follows that

$\frac{v_{\mathrm{r}}}{v}\,\mathbf{e}_{3}\cdot\mathbf{e}_{z}-\frac{v_{z}v_{\mathrm{t}}}{v^{2}}=\sin\varphi\sin
I,$ (47)

with $\varphi$ given in Fig. (12) and $v_{z}$ by equation (30c). Overall, we
finally have

$\displaystyle\langle\Delta L_{z}\rangle$
$\displaystyle=\frac{L_{z}}{v}\langle\Delta v_{\parallel}\rangle,$ (48a)
$\displaystyle\langle(\Delta L_{z}^{2}\rangle$
$\displaystyle=\bigg{(}\frac{L_{z}}{L}\bigg{)}^{2}\bigg{(}\frac{L^{2}}{v^{2}}\langle(\Delta
v_{\parallel})^{2}\rangle\\!+\\!\frac{1}{2}\frac{r^{2}v_{\mathrm{r}}^{2}}{v^{2}}\langle(\Delta
v_{\perp})^{2}\rangle\bigg{)}$ (48b)
$\displaystyle+\frac{r^{2}\sin^{2}\varphi}{2}\bigg{(}1\\!-\\!\frac{L_{z}^{2}}{L^{2}}\bigg{)}\langle(\Delta
v_{\perp})^{2}\rangle,$ $\displaystyle\langle\Delta E\Delta L_{z}\rangle$
$\displaystyle=L_{z}\langle(\Delta v_{\parallel})^{2}\rangle,$ (48c)
$\displaystyle\langle\Delta L\Delta L_{z}\rangle$
$\displaystyle=\frac{L_{z}}{L}\bigg{(}\frac{L^{2}}{v^{2}}\langle(\Delta
v_{\parallel})^{2}\rangle+\frac{1}{2}\frac{r^{2}v_{\mathrm{r}}^{2}}{v^{2}}\langle(\Delta
v_{\perp})^{2}\rangle\bigg{)}.$ (48d)

Equations (44) and (48) 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 w.r.t.
$\varphi$ must be performed with care since $L_{z}^{\prime}$ depends on
$\varphi$ (equation 42). Fortunately, taking advantage of the particular
structure of the diffusion coefficients in $E$, $L$ and $L_{z}$ (see equations
44 and 48), we only need to compute the following four integrals over
$\varphi$

$\displaystyle\int\frac{{\mathrm{d}}\varphi}{2\pi}\langle\Delta
v_{\parallel}\rangle,\quad\int\frac{{\mathrm{d}}\varphi}{2\pi}\langle(\Delta
v_{\parallel})^{2}\rangle,$
$\displaystyle\int\frac{{\mathrm{d}}\varphi}{2\pi}\langle(\Delta
v_{\perp})^{2}\rangle,\quad\int\frac{{\mathrm{d}}\varphi}{2\pi}\cos
2\varphi\,\langle(\Delta v_{\perp})^{2}\rangle,$ (49)

where we used ${\sin^{2}\varphi\\!=\\!\frac{1}{2}(1\\!-\\!\cos 2\varphi)}$ for
convenience.

In equation (3), the only dependence w.r.t. $\varphi$ is in the dependence
w.r.t. $L_{z}^{\prime}$ in $F_{\mathrm{rot}}$. In the particular case of the
LBD (equation 5), we can perform explicitly the average over $\varphi$.
Computing the integrals (49) only requires the evaluation of the two non-
trivial integrals

$\displaystyle A_{1}$
$\displaystyle=\int_{0}^{2\pi}\frac{{\mathrm{d}}\varphi}{2\pi}\,\mathrm{sgn}(L_{z}^{\prime}/L^{\prime}),$
(50a) $\displaystyle A_{2}$
$\displaystyle=\int_{0}^{2\pi}\frac{{\mathrm{d}}\varphi}{2\pi}\,\cos
2\varphi\,\mathrm{sgn}(L_{z}^{\prime}/L^{\prime}),$ (50b)

with $L_{z}^{\prime}$ following from equation (42). Introducing
${\zeta\\!=\\!-(v_{1}^{\prime}v_{\mathrm{t}}/v\\!-\\!v_{3}^{\prime}v_{\mathrm{r}}/v)\cos
I}$ and ${\mu\\!=\\!\zeta/(v_{2}^{\prime}\sin I)}$, these two integrals can be
explicitly computed. Their expressions are gathered in Table 1.

$v_{2}^{\prime}$ | $\mu$ | $A_{1}$ | $A_{2}$
---|---|---|---
0 | $\mathrm{undefined}$ | $-\mathrm{sgn}(r\zeta)$ | 0
$\neq 0$ | $>1$ | $-1$ | 0
$\neq 0$ | $[-1,1]$ | $-2\sin^{-1}\mu/\pi$ | $-2\mu\sqrt{1-\mu^{2}}/\pi$
$\neq 0$ | $<-1$ | $1$ | 0
Table 1: Explicit values of $A_{1}$ and $A_{2}$ (equation 50a) for the LBD DF
(equation 5).

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

$D_{X}=\frac{\Omega_{r}}{\pi}\int_{r_{\mathrm{p}}}^{r_{\mathrm{a}}}\frac{{\mathrm{d}}r}{|v_{\mathrm{r}}|}\langle\Delta
X\rangle_{\varphi}(r),$ (51)

where $\langle\cdot\rangle_{\varphi}$ stands for the angle average. As
described in appendix F2 of T22, equation (51) can be efficiently performed by
changing variable from the radius to an effective anomaly. In particular, this
removes the integrable singularity $1/|v_{\mathrm{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 $\bm{x}$, the FP equation generically reads (Risken,
1996)

$\displaystyle\frac{\partial F(\bm{x},t)}{\partial t}$
$\displaystyle=-\frac{\partial}{\partial\bm{x}}\cdot\bm{\mathcal{F}}_{\bm{x}}(\bm{x})$
(52)
$\displaystyle=-\frac{\partial}{\partial\bm{x}}\cdot\bigg{[}\bm{D}_{\bm{x}}(\bm{x})\,F(\bm{x})-\frac{1}{2}\frac{\partial}{\partial\bm{x}}\\!\cdot\\!\bigg{(}\bm{D}_{\bm{x}\bm{x}}(\bm{x})\,F(\bm{x})\bigg{)}\bigg{]},$

with the diffusion coefficients

$\bm{D}_{\bm{x}}(\bm{x})\\!=\\!\begin{bmatrix}D_{x_{1}}\\\ D_{x_{2}}\\\
D_{x_{3}}\end{bmatrix};\quad\bm{D}_{\bm{x}\bm{x}}(\bm{x})\\!=\\!\begin{bmatrix}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{bmatrix},$
(53)

and the DF in $\bm{x}$ space, ${F(\bm{x})}$. In our case, depending on the
context, ${\bm{x}\\!=\\!(x_{1},x_{2},x_{3})}$ may either stand for
$(E,L,L_{z})$, $(J_{\mathrm{r}},L,L_{z})$ or $(J_{\mathrm{r}},L,\cos I)$.

Following Risken (1996) and Bar-Or & Alexander (2016), one can easily rewrite
the FP equation within some new coordinates ${\bm{x}^{\prime}(\bm{x})}$. The
new diffusion coefficients read

$\displaystyle D_{x_{l}^{\prime}}$ $\displaystyle=\sum_{k}\frac{\partial
x^{\prime}_{l}}{\partial
x_{k}}D_{x_{k}}+\frac{1}{2}\sum_{k,r}\frac{\partial^{2}x^{\prime}_{l}}{\partial
x_{r}\partial x_{k}}D_{x_{r}x_{k}},$ (54a) $\displaystyle
D_{x_{l}^{\prime}x_{m}^{\prime}}$ $\displaystyle=\sum_{k,r}\frac{\partial
x^{\prime}_{l}}{\partial x_{r}}\frac{\partial x^{\prime}_{m}}{\partial
x_{k}}D_{x_{r}x_{k}}.$ (54b)

These coefficients source a FP equation in $\bm{x}^{\prime}$-space for the DF,
${F^{\prime}(\bm{x}^{\prime})}$, reading

$F^{\prime}(\bm{x}^{\prime})=\bigg{|}\frac{\partial\bm{x}}{\partial\bm{x}^{\prime}}\bigg{|}\,F(\bm{x}),$
(55)

with $|\partial\bm{x}/\partial\bm{x}^{\prime}|$ the inverse Jacobian of the
coordinate transform.

###  D.2 From $(E,L,L_{z})$ to $(J_{\mathrm{r}},L,L_{z})$

Following section 2, we have at our disposal diffusion coefficients within the
coordinates ${(E,L,L_{z})}$. Owing to equations (54), the diffusion
coefficients in the action space $(J_{\mathrm{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

$\displaystyle\bm{D}_{1}(\bm{J})\\!=\\!\left[\begin{array}[]{l}D_{J_{\mathrm{r}}}\\\
D_{L}\\\
D_{L_{z}}\end{array}\right];\quad\bm{D}_{2}(\bm{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_{LL}\\!\\!\\!&\\!\\!\\!D_{LL_{z}}\\\
D_{J_{\mathrm{r}}L_{z}}\\!\\!\\!&\\!\\!\\!D_{LL_{z}}\\!\\!\\!&\\!\\!\\!D_{L_{z}L_{z}}\end{array}\right]\\!,\\!$
(62)

which source the FP evolution of $F_{\mathrm{rot}}(J_{\mathrm{r}},L,L_{z})$,
the cluster’s DF in $(J_{\mathrm{r}},L,L_{z})$.

###  D.3 From $(J_{\mathrm{r}},L,L_{z})$ to $(J_{\mathrm{r}},L,\cos I)$

Let us now change of coordinates from
${\bm{J}\\!=\\!(J_{\mathrm{r}},L,L_{z})}$ to
${\bm{J}_{\mathrm{c}}\\!=\\!(J_{\mathrm{r}},L,\cos I\\!=\\!L_{z}/L)}$.
Applying equations (54), we obtain the new diffusion coefficients

$\displaystyle D_{\cos\\!I}$ $\displaystyle=-\frac{\cos
I\,D_{L}}{L}\\!+\\!\frac{D_{L_{z}}}{L}\\!+\\!\frac{\cos
I\,D_{LL}}{L^{2}}\\!-\\!\frac{D_{LL_{z}}}{L^{2}},$ (63a) $\displaystyle
D_{J_{\mathrm{r}}\cos\\!I}$ $\displaystyle=-\frac{\cos
I\,D_{J_{\mathrm{r}}L}}{L}\\!+\\!\frac{D_{J_{\mathrm{r}}L_{z}}}{L}=0,$ (63b)
$\displaystyle D_{L\cos\\!I}$ $\displaystyle=-\frac{\cos
I\,D_{LL}}{L}\\!+\\!\frac{D_{LL_{z}}}{L}=0,$ (63c) $\displaystyle
D_{\cos\\!I\\!\cos\\!I}$
$\displaystyle=\frac{\cos^{2}I\,D_{LL}}{L^{2}}\\!-\\!\frac{2\cos
I\,D_{LL_{z}}}{L^{2}}\\!+\\!\frac{D_{L_{z}L_{z}}}{L^{2}},$ (63d)

while the other coefficients stay unchanged. Importantly, we point out that
${D_{J_{\mathrm{r}}\cos I}\\!=\\!D_{L\cos I}\\!=\\!0}$. This comes from the
relations ${D_{EL_{z}}\\!=\\!\cos I\,D_{EL}}$ and ${D_{LL_{z}}\\!=\\!\cos
I\,D_{LL}}$, which are inferred from equations (44) and (48). Hence, we are
left with the diffusion coefficients

$\displaystyle\bm{D}_{1}(\bm{J}_{\mathrm{c}})\\!=\\!\left[\begin{array}[]{l}D_{J_{\mathrm{r}}}\\\
D_{L}\\\
D_{\cos\\!I}\end{array}\right];\quad\bm{D}_{2}(\bm{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_{LL}\\!\\!\\!&\\!\\!\\!0\\\
0\\!\\!\\!&\\!\\!\\!0\\!\\!\\!&\\!\\!\\!D_{\cos\\!I\\!\cos\\!I}\end{array}\right]\\!.\\!$
(70)

Overall, they drive the FP evolution of
${F(\bm{J}_{\mathrm{c}})\\!=\\!LF_{\mathrm{rot}}(\bm{J})}$, the cluster’s DF
in ${(J_{\mathrm{r}},L,\cos I)}$.

##  Appendix E $N$-body simulations

The simulations presented throughout the main text are performed using the
direct $N$-body code NBODY6++GPU (Wang et al., 2015), version 4.1. We use the
exact same run parameters as in appendix G of T22, and draw 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 2 for all
our binning parameters.

$q$ | 1 | 0 | -6
---|---|---|---
$\small N_{\mathrm{run}}$ | 50 | 50 | 50
$\small t_{\mathrm{last}}\,[{\mathrm{HU}}]$ | 1000 | 1000 | 1000
$\small(J_{\mathrm{r}}^{\min},J_{\mathrm{r}}^{\max})$ | (0, 0.55) | (0, 0.55) | (0, 0.5)
$(\small L^{\min},L^{\max})$ | (0,1.05) | (0,1.05) | (0,1.1)
$\small(N_{J_{\mathrm{r}}},N_{L})$ | (20,20) | (20,20) | (30,40)
$\small(N_{J_{\mathrm{r}}},N_{\cos I})$ | (20,20) | (20,20) | (20,20)
$\small(N_{L},N_{\cos I})$ | (20,20) | (20,20) | (20,20)
Table 2: Detailed parameters for the measurements in $N$-body simulations,
following the same notation as in appendix G1 of T22. To measure relaxation
rate in Fig. 7, we bin the ${(J_{\mathrm{r}},L)}$ domain in
${N_{J_{\mathrm{r}}}\\!\times\\!N_{L}}$ uniform bins within the region
${J_{\mathrm{r}}^{\min}\\!\leq\\!J_{\mathrm{r}}\\!\leq\\!J_{\mathrm{r}}^{\max}}$
(similarly for $L$). We use a similar approach for ${(J_{\mathrm{r}},\cos I)}$
and ${(L,\cos I)}$. All quantities are in physical units
${G\\!=\\!M\\!=\\!b\\!=\\!1}$, if not stated otherwise.

Each $N$-body realisation was composed of ${N\\!=\\!10^{5}}$ (resp.
${N\\!=\\!10^{4}}$) stars and integrated up to
${t_{\max}\\!=\\!1000\,{\mathrm{HU}}}$ (resp.
${t_{\max}\\!=\\!4000\,{\mathrm{HU}}}$) with a dump every ${1{\mathrm{HU}}}$.
On a 40-core CPU node with a single V100 GPU, one simulation typically
required ${\sim\\!24}$ h (resp. ${\sim\\!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_{\mathrm{r}}$; (ii) the determination of the
bins’ size in action space; (iii) the estimation of the relaxation rate
$\partial F/\partial 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_{\mathrm{r}}\\!=\\!0}$ axis. For these three
difficulties, we use the same approach as in T22.

##  Appendix F Sphericality

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

$\bm{\mathcal{I}}=\sum_{k=1}^{N}\rho_{k}^{2}\bigg{(}\bm{r}_{k}^{\mathrm{T}}\bm{r}_{k}\mathbf{I}-\bm{r}_{k}\bm{r}_{k}^{\mathrm{T}}\bigg{)}\bigg{/}\sum_{k=1}^{N}\rho_{k}^{2},$
(71)

with $\bm{r}_{k}$ the location of the $k$-th particle, $\rho_{k}$ its local
density (see Casertano & Hut, 1985) and $\mathbf{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 $\bm{\mathcal{I}}$ is obviously symmetric. It is also semi-definite
positive since for any $\mathbf{y}\in\mathbb{R}^{3}$, we have

$\displaystyle\mathbf{y}^{\mathrm{T}}\bm{\mathcal{I}}\,\mathbf{y}$
$\displaystyle=\sum_{k}\rho_{k}^{2}\,\bigg{(}\bm{r}_{k}^{\mathrm{T}}\bm{r}_{k}\mathbf{y}^{\mathrm{T}}\mathbf{y}-\mathbf{y}^{\mathrm{T}}\bm{r}_{k}\bm{r}_{k}^{\mathrm{T}}\mathbf{y}\bigg{)}\bigg{/}\sum_{k}\rho_{k}^{2}$
$\displaystyle=\sum_{k}\rho_{k}^{2}\,\bigg{(}|\bm{r}_{k}|^{2}\,|\mathbf{y}_{2}|^{2}\\!-\\!|\bm{r}_{k}\\!\cdot\\!\mathbf{y}|^{2}\bigg{)}\bigg{/}\sum_{k}\rho_{k}^{2}\geq
0,$ (72)

following Cauchy–Schwarz’s inequality. As a result, $\bm{\mathcal{I}}$ has
three positive eigenvalues, ${\\{\lambda_{i}\\}_{i}}$, which encapsulate the
cluster’s sphericality. Indeed, spherically symmetric clusters have all their
eigenvalues equal.

We generically define the cluster’s sphericality via
${h\\!=\\!\lambda_{\min}/\lambda_{\max}}$, which we estimate from $N$-body
simulations. To reduce shot noise in that measurement, we average $h$ over
realisations as follows. First, we compute the elementary symmetric
polynomials ${\alpha\\!=\\!\lambda_{1}\\!+\\!\lambda_{2}\\!+\\!\lambda_{3}}$,
${\beta\\!=\\!\lambda_{1}\lambda_{2}\\!+\\!\lambda_{1}\lambda_{3}\\!+\\!\lambda_{2}\lambda_{3}}$
and ${\gamma\\!=\\!\lambda_{1}\lambda_{2}\lambda_{3}}$ for every cluster and
every timestep. We then average the values of ${(\alpha,\beta,\gamma)}$ over
all realisations. From these values, we estimate the eigenvalues $\lambda_{i}$
as the three (positive) roots of the polynomial
${\lambda^{3}\\!-\\!\langle\alpha\rangle\lambda^{2}\\!+\\!\langle\beta\rangle\lambda\\!-\\!\langle\gamma\rangle}$.333The
positivity of $\alpha,\beta,\gamma$ – and hence of their averages – ensures
the positivity of these roots.

In Fig. 2, we illustrate the evolution of the sphericality, $h$, for various
amounts of anisotropy and rotation. We typically find ${h\\!\simeq\\!0.996}$,
i.e. clusters remain reasonably spherically symmetric throughout their
evolution. As as visual check of this conclusion, we represent in Fig. 13 the
late time stellar distribution of a rotating isotropic globular cluster with
${N\\!=\\!10^{4}}$. Even close to core collapse, the cluster remains
spherically symmetric.

Figure 13: Snapshots of the distribution of stars in a rotating
(${\alpha\\!=\\!0.25}$) isotropic (${q\\!=\\!0}$) Plummer cluster with
${N\\!=\\!10^{4}}$ stars, projected on the rotation plane. From left to right,
this corresponds to ${t\\!=\\!0,1000,2000,3000\,{\mathrm{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

$F_{\mathrm{rot}}(E,L_{z})\propto(\mathrm{e}^{-\beta E}-1)\
\mathrm{e}^{-\beta\Omega_{0}L_{z}},$ (73)

where
${F_{\mathrm{rot}}\\!=\\!F_{\mathrm{rot}}(\bm{r},\bm{v})\\!=\\!F_{\mathrm{rot}}(E,L_{z})}$
is the DF in ${(\bm{r},\bm{v})}$. In this DF,
${\beta\\!=\\!1/\sigma_{\mathrm{c}}^{2}}$ is an inverse temperature with
$\sigma_{\mathrm{c}}$ the central velocity dispersion, while
${\omega_{0}\\!=\\!\sqrt{9/(4\pi Gn_{\mathrm{c}})}\,\Omega_{0}}$ is a rotation
parameter with $n_{\mathrm{c}}$ the central density and $\Omega_{0}$ the
angular velocity in the cluster’s centre (Lagoute & Longaretti, 1996). In
particular, $\beta$ is related to the so-called King parameter (King, 1966),
defined by ${W_{0}\\!=\\!-\beta(\psi-\psi_{\mathrm{t}})}$ where
$\psi_{\mathrm{t}}$ is the potential at the cluster’s edge. The rotating King
models from Einsel & Spurzem (1999) can therefore be parametrised by
$(W_{0},\omega_{0})$. In practice, Einsel & Spurzem (1999) integrated these
models up to core collapse using a FP scheme in ${(E,L_{z})}$ for different
rotation parameters, $\omega_{0}$. 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 $\omega_{0}$, 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_{\mathrm{r}},L,L_{z})}$
(section D.2) or in ${(J_{\mathrm{r}},L,\cos I)}$ (section D.3). In this
appendix, we compute the 2D FP equations in ${(J_{\mathrm{r}},L)}$,
${(J_{\mathrm{r}},\cos I)}$ and ${(L,\cos I)}$ that respectively drive the
evolution of the DFs

$\displaystyle F(J_{\mathrm{r}},L)$
$\displaystyle=\int_{-L}^{L}{\mathrm{d}}L_{z}\,F_{\mathrm{rot}}(J_{\mathrm{r}},L,L_{z}),$
(74a) $\displaystyle F(J_{\mathrm{r}},\cos I)$
$\displaystyle=\int_{0}^{+\infty}{\mathrm{d}}L\,F(J_{\mathrm{r}},L,\cos I),$
(74b) $\displaystyle F(L,\cos I)$
$\displaystyle=\int_{0}^{+\infty}{\mathrm{d}}J_{\mathrm{r}}\,F(J_{\mathrm{r}},L,\cos
I).$ (74c)

###  H.1 Equation in $(J_{\mathrm{r}},L)$

We start from equation (1), which describes diffusion in
${(J_{\mathrm{r}},L,L_{z})}$, and write the flux divergence as

$\displaystyle\frac{\partial F_{\mathrm{rot}}}{\partial t}$
$\displaystyle=-\bigg{(}\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}}\bigg{)}.$ (75)

Integrating over $L_{z}$ yields

$\displaystyle\\!\\!\int_{-L}^{L}\\!\\!\\!\\!{\mathrm{d}}L_{z}\,\frac{\partial\mathcal{F}_{J_{\mathrm{r}}}}{\partial
J_{\mathrm{r}}}$ $\displaystyle=\frac{\partial}{\partial
J_{\mathrm{r}}}\\!\\!\int_{-L}^{L}\\!\\!\\!\\!{\mathrm{d}}L_{z}\,\mathcal{F}_{J_{\mathrm{r}}},$
(76a)
$\displaystyle\\!\\!\int_{-L}^{L}\\!\\!\\!\\!{\mathrm{d}}L_{z}\,\frac{\partial\mathcal{F}_{L}}{\partial
L}$
$\displaystyle=-\mathcal{F}_{L}(L_{z}\\!=\\!L)\\!-\\!\mathcal{F}_{L}(L_{z}\\!=\\!-L)+\frac{\partial}{\partial
L}\\!\\!\int_{-L}^{L}\\!\\!\\!\\!{\mathrm{d}}L_{z}\,\mathcal{F}_{L},$ (76b)
$\displaystyle\\!\\!\int_{-L}^{L}\\!\\!\\!\\!{\mathrm{d}}L_{z}\,\frac{\partial\mathcal{F}_{L_{z}}}{\partial
L_{z}}$
$\displaystyle=\mathcal{F}_{L_{z}}(L_{z}\\!=\\!L)-\mathcal{F}_{L_{z}}(L_{z}\\!=\\!-L).$
(76c)

In addition, because the flux cannot exit action space, we also have
${\mathcal{F}_{L_{z}}(L_{z}\\!=\\!\pm
L)\\!=\\!\pm\mathcal{F}_{L}(L_{z}\\!=\\!\pm L)}$. As a result, all the non-
integral terms in equations (76) cancel one another. We are then left with the
2D equation

$\frac{\partial F(J_{\mathrm{r}},L,t)}{\partial
t}=-\frac{\partial}{\partial(J_{\mathrm{r}},L)}\cdot\\!\\!\int_{-L}^{L}\\!\\!{\mathrm{d}}L_{z}\,\begin{bmatrix}\mathcal{F}_{J_{\mathrm{r}}}(J_{\mathrm{r}},L,L_{z})\\\\[4.30554pt]
\mathcal{F}_{L}(J_{\mathrm{r}},L,L_{z})\end{bmatrix}.$ (77)

It is this rewriting that is used in section 4.2. To obtain Fig. 8, the
$L_{z}$-integrals in equations (76) are sampled with 50 nodes. The
${(w,\vartheta,\phi)}$-integrals of equations (3) are sampled with 100
$w$-nodes, 100 $\vartheta$-nodes and 200 $\phi$-nodes.

###  H.2 Equation in $(J_{\mathrm{r}},\cos I)$

We start from the diffusion coefficients in equation (70), which describe
diffusion in ${(J_{\mathrm{r}},L,\cos I)}$, for the DF,
${F\\!=\\!LF_{\mathrm{rot}}}$. We write the flux divergence as

$\displaystyle\frac{\partial F}{\partial t}$
$\displaystyle=-\bigg{(}\frac{\partial\mathcal{F}_{J_{\mathrm{r}}}}{\partial
J_{\mathrm{r}}}+\frac{\partial\mathcal{F}_{L}}{\partial
L}+\frac{\partial\mathcal{F}_{\cos\\!I}}{\partial\cos\\!I}\bigg{)}.$ (78)

Integrating over $L$ yields

$\displaystyle\\!\\!\int_{0}^{\infty}\\!\\!\\!\\!{\mathrm{d}}L\,\frac{\partial\mathcal{F}_{J_{\mathrm{r}}}}{\partial
J_{\mathrm{r}}}$ $\displaystyle=\frac{\partial}{\partial
J_{\mathrm{r}}}\\!\\!\int_{0}^{+\infty}\\!\\!\\!\\!{\mathrm{d}}L\,\mathcal{F}_{J_{\mathrm{r}}},$
(79a)
$\displaystyle\\!\\!\int_{0}^{\infty}\\!\\!\\!\\!{\mathrm{d}}L\,\frac{\partial\mathcal{F}_{L}}{\partial
L}$
$\displaystyle=\mathcal{F}_{L}(L\\!=\\!+\infty)-\mathcal{F}_{L}(L\\!=\\!0),$
(79b)
$\displaystyle\\!\\!\int_{0}^{\infty}\\!\\!\\!\\!{\mathrm{d}}L\,\frac{\partial\mathcal{F}_{\cos\\!I}}{\partial\cos\\!I}$
$\displaystyle=\frac{\partial}{\partial\cos\\!I}\\!\\!\int_{0}^{+\infty}\\!\\!\\!\\!{\mathrm{d}}L\,\mathcal{F}_{\cos\\!I}.$
(79c)

Because the flux cannot leave action space, equation (79b) vanishes. Then,
equation (78) readily reduces to a 2D FP equation in $(J_{\mathrm{r}},\cos
I)$, which we use in section 5.2. To obtain Fig. 10, the $L$-integrals in
equations (79) are sampled with 50 nodes for
${0\\!\leq\\!L\\!\leq\\!3\,L_{0}}$, where ${L_{0}\\!=\\!\sqrt{GMb}}$ is the
typical action. The ${(w,\vartheta,\phi)}$-integrals of equations (3) are
sampled with 100 $w$-nodes, 100 $\vartheta$-nodes and 800 $\phi$-nodes.

###  H.3 Equation in $(L,\cos I)$

As in appendix H.2, we start from the diffusion coefficients in equation (70).
Integrating equation (78) over $J_{\mathrm{r}}$ yields

$\displaystyle\\!\\!\int_{0}^{\infty}\\!\\!\\!\\!{\mathrm{d}}J_{\mathrm{r}}\,\frac{\partial\mathcal{F}_{J_{\mathrm{r}}}}{\partial
J_{\mathrm{r}}}$
$\displaystyle=\mathcal{F}_{J_{\mathrm{r}}}(J_{\mathrm{r}}\\!=\\!+\infty)-\mathcal{F}_{J_{\mathrm{r}}}(J_{\mathrm{r}}\\!=\\!0),$
(80a)
$\displaystyle\\!\\!\int_{0}^{\infty}\\!\\!\\!\\!{\mathrm{d}}J_{\mathrm{r}}\,\frac{\partial\mathcal{F}_{L}}{\partial
L}$ $\displaystyle=\frac{\partial}{\partial
L}\\!\\!\int_{0}^{+\infty}\\!\\!\\!\\!{\mathrm{d}}J_{\mathrm{r}}\,\mathcal{F}_{L}\,,$
(80b)
$\displaystyle\\!\\!\int_{0}^{\infty}\\!\\!\\!\\!{\mathrm{d}}J_{\mathrm{r}}\,\frac{\partial\mathcal{F}_{\cos\\!I}}{\partial\cos\\!I}$
$\displaystyle=\frac{\partial}{\partial\cos\\!I}\\!\\!\int_{0}^{+\infty}\\!\\!\\!\\!{\mathrm{d}}J_{\mathrm{r}}\,\mathcal{F}_{\cos\\!I}.$
(80c)

Once again, because the flux cannot leave action space, equation (80a)
vanishes. Then, equation (78) reduces to a 2D FP equation in $(L,\cos I)$,
which we use in appendix I. To obtain Fig. 15, the $J_{\mathrm{r}}$-integrals
in equations (80) are sampled with 50 nodes for
${0\\!\leq\\!J_{\mathrm{r}}\\!\leq\\!10\,L_{0}}$. The
${(w,\vartheta,\phi)}$-integrals of equations (3) are sampled with 100
$w$-nodes, 100 $\vartheta$-nodes and 1600 $\phi$-nodes.

##  Appendix I Relaxation in $(L,\cos I)$

In this appendix, we follow the same approach as in section 5, and investigate
relaxation in ${(L,\cos I)}$. We refer to Appendix H.3 for the derivation of
the relevant 2D FP equation. In Fig. 14, we illustrate the $N$-body
measurements while Fig. 15 presents the associated NR predictions.

Figure 14: Same as in Fig. 9 but in $(L,\cos I)$. Diffusion reshuffles orbital
inclinations toward a more affine distribution in $\cos I$.  Figure 15: Same
as Fig. 10 but in $(L,\cos I)$. The NR prediction fails to recover in detail
the diffusion structures observed numerically in Fig. 14.

From these two figures, we can globally make the same observations as in
section 5, which considered diffusion in ${(J_{\mathrm{r}},\cos I)}$. Indeed,
in rotating clusters, diffusion reshuffles orbits towards smoother
distributions of inclinations. Nevertheless, the NR prediction fails at
predicting the exact structures measured in $N$-body simulations. In
particular, the NR theory predicts a relaxation rate varying weakly with $\cos
I$ (for a given sign of $\cos I$). This differs from the $N$-body simulations,
where the relaxation rate decreases away from ${\cos I\\!=\\!0}$.

##  Appendix J Impact of discontinuities

In this appendix, we investigate in more detail the impact of the
discontinuity at ${\cos I\\!=\\!0}$, introduced by the LBD (equation 5). To do
so, we follow the same approach as in section 2.3 of Rozier et al. (2019) and
consider rotating DFs of the form

$F_{\mathrm{rot}}(J_{\mathrm{r}},L,L_{z})=F_{\mathrm{tot}}(J_{\mathrm{r}},L)(1+\alpha
g[L_{z}/L]),$ (81)

where $g[\cos I]$ is an odd function with ${g(1)\\!=\\!1}$. The LBD
corresponds to ${g\\!=\\!\mathrm{sgn}}$. To approximate smoothly the LBD, we
consider the sequence of functions

$g_{a}(x)=\mathrm{erf}(ax)/\mathrm{erf}(a).$ (82)

As illustrated in Fig. 16, this ensures that ${g_{0}(x)\\!=\\!x}$ and
${g_{\infty}(x)\\!=\\!\mathrm{sgn}(x)}$.

Figure 16: Family of functions $g_{a}(\cos I)$ (equation 82) for various $a$.
For $a\\!\rightarrow\\!0$, the function approaches identity, while for
${a\\!\rightarrow\\!\infty}$, it approaches the sign function. Varying $a$
allows us to investigate the impact of the LBD discontinuity in equation (5).

To probe the possible presence of discontinuities and singularities, let us
first compute the $\cos I$ component of the 3D flux in
${(J_{\mathrm{r}},L,\cos I)}$, as introduced in section D.3. The dependence of
this flux w.r.t. $a$ is illustrated in Fig. 17.

Figure 17:  Diffusion flux along $\cos I$ of the 3D FP equation in
${(J_{\mathrm{r}},L,\cos I)}$, for ${J_{\mathrm{r}}\\!=\\!0.1}$ and ${\cos
I\\!=\\!0.2}$, as a function of $L$ and for various smooth functions $g_{a}$
(equation 82). Here, we consider an isotropic cluster with rotation parameter
${\alpha\\!=\\!0.25}$. For any smooth $g_{a}$, the flux diverges like ${1/L}$
for ${L\\!\to\\!0}$. As ${a\\!\to\\!\infty}$, i.e. as $g_{a}$ tends to the
$\mathrm{sgn}$ function, the flux converges to the LBD flux pointwise, which
exhibits no divergence.

In this figure, for any smooth $g_{a}$, we observe a $1/L$ divergence of the
flux as ${L\\!\to\\!0}$. In a nutshell, when working with ${\cos I}$, we
suffer from a coordinate singularity, and the NR prediction cannot be applied
to a DF with a smooth rotation function. Yet, when
${a\\!\rightarrow\\!\infty}$, the flux converges towards the LBD flux
pointwise, and this flux does not diverge for ${L\\!\to\\!0}$. Phrased
differently, in the particular case of the LBD sign function, we can make a
meaningful and well-posed NR prediction for the diffusion in ${\cos I}$. In
practice, this vanishing of the divergence stems from the cancellation of the
derivative of ${\mathrm{sgn}(\cos I)}$ everywhere (except for ${\cos
I\\!=\\!0}$). Such a property is not the norm for smooth arbitrary rotation
functions, ${g(\cos I)}$. In that case, the $1/L^{2}$ singularities visible in
equations (63) do not combine into an integrable quantity.

To further stress that this divergence originate from coordinate
singularities, let us now produce NR predictions by describing relaxation in
${(J_{\mathrm{r}},L,L_{z})}$ rather than in ${(J_{\mathrm{r}},L,\cos I)}$. To
do so, we define the DF in $L_{z}$ as

$F(L_{z})=\int_{0}^{+\infty}\\!\\!\\!\\!{\mathrm{d}}J_{\mathrm{r}}\,\\!\\!\int_{|L_{z}|}^{+\infty}\\!\\!\\!\\!{\mathrm{d}}L\,F_{\mathrm{rot}}(J_{\mathrm{r}},L,L_{z}).$
(83)

Integrating the 3D FP equation (1) over $J_{\mathrm{r}}$ and $L$ yields

$\frac{\partial F(L_{z})}{\partial t}=-\frac{\partial}{\partial
L_{z}}\big{(}\overline{D}_{L_{z}}F[L_{z}]\big{)}+\frac{1}{2}\frac{\partial^{2}}{\partial
L_{z}^{2}}\big{(}\overline{D}_{L_{z}L_{z}}F[L_{z}]\big{)}.$ (84)

Here, the 1D diffusion coefficients in $L_{z}$ are given by

$\displaystyle\overline{D}_{L_{z}}$
$\displaystyle=\\!\\!\int_{0}^{+\infty}\\!\\!\\!\\!{\mathrm{d}}J_{\mathrm{r}}\,\\!\\!\int_{|L_{z}|}^{+\infty}\\!\\!\\!\\!{\mathrm{d}}L\,D_{L_{z}}F_{\mathrm{rot}}(J_{\mathrm{r}},L\,|\,L_{z}),$
(85a) $\displaystyle\overline{D}_{L_{z}L_{z}}$
$\displaystyle=\\!\\!\int_{0}^{+\infty}{\mathrm{d}}J_{\mathrm{r}}\,\\!\\!\int_{|L_{z}|}^{+\infty}\\!\\!\\!\\!{\mathrm{d}}L\,D_{L_{z}L_{z}}F_{\mathrm{rot}}(J_{\mathrm{r}},L\,|\,L_{z}),$
(85b)

with
${F_{\mathrm{rot}}(J_{\mathrm{r}},L\,|\,L_{z})\\!=\\!F_{\mathrm{rot}}(J_{\mathrm{r}},L,L_{z})/F(L_{z})}$
standing for the DF in $(J_{\mathrm{r}},L)$ given $L_{z}$, and normalised to
unity. Importantly, for ${L_{z}\\!\neq\\!0}$, the coefficients
$\overline{D}_{L_{z}}$ and $\overline{D}_{L_{z}L_{z}}$ are well defined
whatever the considered rotating DF. In addition, they both converge to some
finite values as ${L_{z}\\!\to\\!0}$, as illustrated in Fig. 18.

Figure 18:  One-dimensional diffusion coefficients $\overline{D}_{L_{z}}$ (top
panel) and $\overline{D}_{L_{z}L_{z}}$ (bottom panel) of equations (85),
computed for the rotating isotropic clusters with ${\alpha\\!=\\!0.25}$, and
with the smoothing parameters ${a\\!=\\!1,5,10,15,\infty}$, as defined in
equation (82). These coefficients are well-defined for all values of $L_{z}$.
In addition, they only (very) weakly depend on $a$, with relative differences
of order 1–5%.

The absence of any divergence here emphasizes that the divergence observed in
Fig. 17 stems from a coordinate singularity associated with ${\cos I}$.

  *[BH]: black hole
  *[NR]: non-resonant relaxation
  *[DF]: distribution function
  *[FP]: Fokker–Planck
  *[DFs]: distribution function
  *[LBD]: Lynden-Bell daemon
  *[VRR]: vector resonant relaxation
  *[RR]: resonant relaxation
  *[BHs]: black hole