The longterm relaxation of rotating spherically symmetric globular clusters is investigated through an extension of the orbitaveraged Chandrasekhar nonresonant formalism. A comparison is made with the longterm evolution of the distribution function in action space measured from averages of sets of Nbody simulations up to core collapse. The impact of rotation on inplane relaxation is found to be weak. In addition we observe a clear match between theoretical predictions and Nbody measurements. For the class of rotating models considered we find no strong gravogyro 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 nonresonant prediction does not reproduce numerical measurements. We argue that this mismatch stems from these orbits coherent interactions which are not captured by the nonresonant formalism that only addresses local deflections.
['diffusion', 'gravitation', 'galaxies: kinematics and dynamics', 'Astrophysics - Astrophysics of Galaxies']
Nonresonant relaxation of rotating globular clusters
{'Non-resonant relaxation of rotating globular clusters': "Kerwann Tep 1 , 2 ⋆ , Jean-Baptiste Fouvry 1 , and Christophe Pichon 1 , 3 , 4  \n- 1 Institut d'Astrophysique de Paris, CNRS and Sorbonne Université, UMR 7095, 98 bis Boulevard Arago, F-75014 Paris, France\n- 2 Department of Physics and Astronomy, The University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA\n- 3 IPhT, DRF-INP, UMR 3680, CEA, L'Orme des Merisiers, Bât 774, 91191 Gif-sur-Yvette, France\n- 4 Korea Institute for Advanced Study, 85 Hoegi-ro, Dongdaemun-gu, Seoul 02455, Republic of Korea  \nSeptember 16, 2024", 'ABSTRACT': "The long-term relaxation of rotating, spherically symmetric 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 reshu ffl ing 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.  \nKey words. Di ff usion - Gravitation - Galaxies: kinematics and dynamics", '1. Introduction': "Rotation is ubiquitous in stellar systems. In e ff ect, it provides a source of free energy, allowing clusters to e ffi ciently reshu ffl e their orbital structure towards more likely configurations. Yet, historically, the study of globular clusters has been mainly focused on isotropic, non-rotating, old globular clusters (Aarseth et al. 1974; Spitzer 1975; Cohn 1979; Trager et al. 1995; Miocchi et al. 2013). The reasons behind such simplifications are two-fold: (i) naturally, it is easier numerically and analytically to neglect the e ff ect of rotation; (ii) spherical isotropic models for example, the King models (King 1966) or the Wilson models (Wilson 1975) - provided a satisfactory zeroth-order description of the main observed dynamical properties of globular clusters (see, e.g., McLaughlin & van der Marel 2005).  \nThe last decade has seen the extraction of new data - for instance, 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.  \nTherefore, a secular theory that describes the evolution of globular clusters and accounts for their rotation is needed to describe their long-term evolution. In particular, the Fokker-Planck theory was used to probe the impact of rotation of the long-term evolution of such clusters, both before (Einsel & Spurzem 1999) and after (Kim et al. 2002) core collapse. Furthermore, observations show that the angular momentum distribution measured  \nin 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).  \nIn isolated systems, a non-zero total angular momentum, that is, the presence of rotation, can have a significant impact on the cluster's long-term evolution, for instance, through the 'gravo-gyro catastrophe' (Hachisu 1979; Ernst et al. 2007). This phenomenon has been observed in a range of rotating systems, including gas cylinders (Inagaki & Hachisu 1978; Hachisu 1979), gaseous discs (Hachisu 1982), flattened (quasi-spherical) star clusters (Akiyama & Sugimoto 1989; Einsel & Spurzem 1999; Ernst et al. 2007), and clusters with embedded black holes (BHs) (Fiestas & Spurzem 2010; Kamlah et al. 2022). Let us however note that many of these studies used a rotating King model to study the impact of rotation (see, e.g., Einsel & Spurzem 1999; Varri & Bertin 2012). Therein, changing the amount of rotation impacts the density profile: this makes comparisons amongst di ff erent models less clear.  \nThe concurrent occurrence of internal rotation and a spectrum of stellar masses can result in the formation of an oblate core of fast rotating heavy masses (Kim et al. 2004; Tiongco et al. 2021). More precisely, the orbital inclinations of the heaviest stars align with respect to one another, inducing a mass segregation in the distribution of orbital inclinations (Szölgyén et al. 2019). The generally agreed explanation for this phenomenon is resonant relaxation and resonant friction (Rauch & Tremaine 1996; Meiron & Kocsis 2019). While this e ff ect concerns globular clusters, nuclear cluster with a dominant massive BH (Szölgyén & Kocsis 2018; Foote et al. 2020; Gruzinov et al. 2020;  \nMagnan et al. 2022; Ginat et al. 2023) also display this spontaneous alignment.  \nThe objective of this paper is to explore the impact of rotation on the long-term evolution of spherically symmetric globular clusters. To achieve this goal, we extend the orbit average analysis of anisotropic Plummer globular clusters from Tep et al. (2022) to the case of rotating clusters. By performing tailored N -body simulations, we quantify the validity of Chandrasekhar's non-resonant relaxation (NR) theory (Chandrasekhar 1943), and assess the importance of the gravo-gyro catastrophe as well as that of collective e ff ects (i.e., the self-amplification of the star's gravitational response, see Heyvaerts 2010).  \nThis 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 di ff usion) and 5 (out-of-plane di ff usion). Finally, we discuss our results in section 6.", '2. Non-resonant relaxation': 'We consider a self-gravitating globular cluster composed of N stars of individual mass m = M / N , with M the cluster\'s total mass. We follow the cluster\'s evolution through the total distribution function (DF) in ( r , v ) space, F rot = F rot( r , v ), with r the position and v the velocity. All normalisations are taken to be the same as in section 2 of Tep et al. (2022), hereafter T22. In addition to quasi-stationarity, we assume the cluster to be spherically symmetric - hence, with planar unperturbed orbits - and in rotation. Since we are interested in orbital distortion, it is convenient to monitor the cluster\'s evolution in action space, via F rot = F rot( J ), where J = ( J r , L , Lz ) are the specific action coordinates. Here, J r is the radial action, L the norm of the angular momentum vector and Lz its projection along the z -axis. We also introduce the orbital inclination through cos I = Lz / L . Because actions are integrals of motion, we use them to label the cluster\'s orbits and track their deformation over time.  \nThe long-term evolution of the cluster\'s DF is described by the Fokker-Planck (FP) equation (see, e.g., section 7.4 of Binney &Tremaine 2008)  \n∂ F rot( J , t ) ∂ t = -∂ ∂ J · F ( J ) (1) = -∂ ∂ J · " D 1( J ) F rot( J ) -1 2 ∂ ∂ J · D 2( J ) F rot( J ) !# ,  \nwith F ( J ) the action space flux. Therein, we find the first-order di ff usion coe ffi cient, D 1, and the di ff usion matrix, 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 di ff usion flux in equation (1).', '2.1. Local velocity deflection coefficients': "Let us consider a test star of mass m and velocity v . As a result of the cluster's finite number of constituents, this test star is subject to perturbations around its mean field trajectory, driving an irreversible di ff usion 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  \ndeflection coe ffi cients generically read  \n⟨ ∆ vi ⟩ = -2 A Z d v ' wi w 3 F rot( r , v ' ) , (2a)  \n⟨ ∆ vi ∆ vj ⟩ = A Z d v ' w 2 δ i j -wiwj w 3 F rot( r , v ' ) , (2b)  \nwhere 1 ≤ i , j ≤ 3 run over the three directions of the coordinate system. We also defined A = 4 π mG 2 ln Λ and the relative velocity w = v -v ' . Here, ln Λ stands for the Coulomb logarithm, set to Λ= 0 . 11 N 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 nonrotating but anisotropic cluster, was pursued in T22.  \nIn 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 di ff usion coe ffi -cients. More precisely, we write equations (2) as  \n          ⟨ ∆ v ∥ ⟩ ⟨ ( ∆ v ∥ ) 2 ⟩ ⟨ ( ∆ v ⊥ ) 2 ⟩           = A Z d w d ϑ d ϕ sin ϑ            -2 cos ϑ w sin 2 ϑ w (1 + cos 2 ϑ )            F rot , (3)  \nwhere ∆ v ∥ (resp. ∆ v ⊥ ) is the local velocity deflection along (resp. perpendicular to) the test star's trajectory and we have used polar coordinates ( w , ϑ, ϕ ) with the z -axis parallel to v (see Fig. A.1). Appendix A.1 also details the arguments at which the DF, F rot( r , v ' ) = F rot( E ' , L ' , L ' z ), must be evaluated.  \nEquation (3) is one of the main results of the present work. In particular, equation (3) possesses a few key advantages compared to the expressions given in equations (3) of T22, which are recovered by the present approach: (i) no integrable singular denominators remain; (ii) no gradients of F rot are required during the computation; (iii) the positivity of the second-order coe ffi -cients 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 di ff usion of orbits, we must consider the e ff ect of NR on orbital invariants, namely the energy E and the angular momenta L and Lz . In appendix B, we expand upon T22 and detail how the local di ff usion coe ffi cients in ( E , L , Lz ) may be computed from the local velocity coe ffi cients given in equation (3).  \nOnce the local di ff usion coe ffi cients 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 twodimensional integral spanning the radial range of the test orbit ( r ∈ [ r p , r a]) and its angular phase ( φ ∈ [0 , 2 π ]) within the orbital plane. This is illustrated in Fig. 1.  \nIn practice, computing the orbit average di ff usion coe ffi -cients involves computing expressions of the form  \nDX = Ω r π Z r a r p d r | v r | Z 2 π 0 d φ 2 π ⟨ ∆ X ⟩ ( r , φ ) , (4)  \nwhere X runs over all the di ff usion coe ffi cients to compute, that is, runs over singlets and pairs from { E , L , Lz } . Here, Ω r is the radial frequency and v r the radial velocity that is evaluated along  \nFig. 2. Sphericity, 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 α = 0 . 1 , 0 . 25 , 0 . 5 (ordered from dark to light colors). In agreement with Rozier et al. (2019), some of these clusters are unstable - namely ( q , α ) = (1 , 0 . 5), and in a smaller fashion ( q , α ) = (1 , 0 . 25) - while the others are linearly stable and remain spherically symmetric.  \n<!-- image -->  \nFig. 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 ∈ [ r p , r a]) and its angular range ( φ ∈ [0 , 2 π ]), weighted by the surface density Ω r / (2 π 2 | v r | ). Here, the frame is that of the right panel of Fig. A.2.  \n<!-- image -->  \nthe mean field orbit. We refer to appendix C for detailed expressions. In particular, in appendix C.1, we show how the angular average can be performed explicitly for some classes of rotating models. And, in appendix C.2, we use the same approach as in appendix F2 of T22 to compute the radial average using a stable numerical integration scheme.  \nThe final stage of the calculation is to convert the orbit-averaged di ff usion coe ffi cients in ( E , L , Lz ) into ones in J = ( J r , L , Lz ). This is rather straightforward, as detailed in appendix D. Lastly, these di ff usion coe ffi cients are those used to evaluate the FP equation (1). Along the same lines, one can similarly describe the dynamics within the coordinates J c = ( J r , L , cos I ).', '3. Long-term relaxation': 'We wish to study the impact of rotation on the long-term evolution of rotating, spherically symmetric 1 , 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.  \nFirst, let us describe the classes of clusters considered. The in-plane distribution of orbits (i.e., when projecting in the ( J r , L ) space) follows the same Plummer DFs (see, e.g., Dejonghe 1987) as in section 3.1 of T22. Velocity anisotropy is encoded with the dimensionless parameter q , with q = 0 corresponding to an isotropic distribution, and q > 0 (resp. q < 0) associated with radially (resp. tangentially) anisotropic velocity distributions.  \nLet 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  \nF rot( J ) = F tot( J r , L ) GLYPH<16> 1 + α sgn[ Lz / L ] GLYPH<17> , (5)  \nwhere F tot( J r , L ) is the non-rotating DF and α a dimensionless parameter between 0 and 1. Physically, this parameter corresponds to converting a fraction α of retrograde orbits ( Lz < 0) into prograde orbits ( Lz > 0). Importantly, since Lz 7→ sgn[ Lz / L ] is an odd function, we stress that F tot and F rot generate the exact same potential.  \nWhen considering orbital inclinations via J c = ( J r , L , cos I ) (see appendix D.3), the LBD yields the DF in J c space 2  \nF ( J c) = L F tot( J r , L ) GLYPH<16> 1 + α sgn[cos I ] GLYPH<17> . (6)  \nIntegrating this equation over cos I , one recovers the reduced DF in ( J r , L ), F ( J r , L ) = 2 LF tot( J r , L ), used to describe non-rotating clusters (see, e.g., Hamilton et al. 2018).', '3.1. Sphericity of the cluster': "The NR theory presented in section 2 assumes that the cluster remains spherically symmetric throughout. However, Rozier et al. (2019) showed that (su ffi ciently) rotating clusters can harbour unstable modes (see, e.g., fig. 8 therein). To avoid this complication, we restricted ourselves to (linearly) stable rotating clusters, that is, clusters that would remain spherically symmetric.  \nTo probe the conservation of spherical symmetry, we investigateed the 'sphericity', h , of clusters, as defined in appendix F in terms of the eigenvalue ratio of the density-weighted inertia tensor. This is illustrated in Fig. 2. Except for the radiallyanisotropic and rapidly-rotating cluster ( q , α ) = (1 , 0 . 5), all clusters remain approximately spherical. This is in agreement with the measurements from Rozier et al. (2019), as this particular cluster falls into the region of linear instability. Therefore, we shall focus our interest on parameters ( q , α ) for which the clusters remain spherically symmetric.", '3.2. The gravo-gyro catastrophe': 'The typical relaxation timescale of a globular cluster can be estimated through the half-mass relaxation time (see, e.g., section 14 of Heggie & Hut 2003, for more details), defined by  \nt rh = 0 . 138 N r 3 / 2 h ( GM ) 1 / 2 ln(0 . 11 N ) , (7)  \nFig. 3. Core radius as a function of time, as measured in N -body simulations, with t rh defined in equation (7). Each case has been averaged over 50 realisations with N = 10 4 . Increasing the rotation strength α slightly reduces the time of core collapse. Nevertheless, the impact of rotation (i.e. the gravo-gyro catastrophe) is surely not as pronounced as what was observed in, e.g., Einsel & Spurzem (1999).  \n<!-- image -->  \nwith r h the half-mass radius. For a cluster with N = 10 5 stars, this yields t rh = 994 HU, where HU stands for the Hénon units (Heggie & Mathieu 1986) in which G = M = R v = 1, with R v the virial radius 3 . In practice, with N = 10 5 running a numerical simulation up to t ∼ t rh took about one day of computation, see appendix E. In Breen et al. (2017), core collapse occurs at t ∼ 17 t rh for an isotropic cluster (see table 1 therein for the dependence with respect to velocity anisotropy). Hence, integrating clusters with N = 10 5 up to core collapse is not reasonably feasible. In this section we therefore scaled down the size of the clusters to N = 10 4 stars. In that case, t rh = 132 HU, and core collapse was numerically reached in about ten hours. Since we focused on measuring the core radius - a very integrated quantity - the quality of the numerical measurements was not too much degraded by this use of a smaller value for N . In practice, measurements were averaged over 50 runs.  \nFigure 3 represents the time evolution of the averaged core radius, R c, defined as (see, e.g., Breen et al. 2017; T22)  \nR 2 c = X i r 2 i ρ 2 i / X i ρ 2 i , (8)  \nwhere ri is the radial position of star i and ρ i an estimator of the density at ri , as defined in NBODY6 ++ GPU (Wang et al. 2015). for the usual sets of anisotropies and rotations ( q = 1 , 0 , -6 and α = 0 , 0 . 1 , 0 . 25 , 0 . 5). Interestingly, for these rotating clusters, the impact of rotation is definitely not as important as the one reported in Einsel & Spurzem (1999) in rotating King models (see figure 2 therein). Appendix G reproduces the main parameters characterising these models. For example, in Einsel & Spurzem (1999), assuming W 0 = 6 . 0, a non-rotating King cluster collapses at t ∼ 12 t rh, while a rotating one, with ω 0 / Ω 0 = 0 . 4, can collapse as early as t ∼ 9 t rh. Here, for Plummer spheres, we do not find any such stark impact of rotation.  \nNonetheless, Figure 3 does not contradict the measurements from Einsel & Spurzem (1999). Indeed, the parameters ω 0 (King) and α (Plummer) do not parametrise rotation in the same way. As Einsel & Spurzem (1999) varies ω 0, the mean density profile gets modified. Indeed, the rotating DF of the King model (equation G.1) cannot be decomposed into a fixed, rotation-free even part, and an odd part in Lz (see, e.g., Dejonghe 1986). As  \nsuch, we argue that the LBD approach ensures a fairer comparison between models to isolate the distinctive impact of rotation.  \nOf course, to better assess the universality of these observations, it would be worthwhile to perform similar experiments with other parametrisations for the rotation (see Appendix K). This is left for future investigations.', '3.3. In-plane vs out-of-plane diffusion': "Using once again N -body simulations with N = 10 4 stars, let us finally investigate the typical features of the in-plane di ff usion (i.e. relaxation in ( J r , L )) and out-of-plane di ff usion (i.e. relaxation in cos I ).  \nFirst, Fig. 4 reports on the time evolution of the DF in L , that is, the in-plane relaxation of the distribution of angular momenta. There, we recover the imprints of core collapse visible through the slow overall contraction of the distribution. But importantly, we clearly note that rotation only (very) weakly impacts this relaxation.  \nWe can now compare this relaxation with the out-of-plane one. Figure 5 reports on the time evolution of the DF in cos I . As could have been expected, the clusters tend towards smoother distribution of inclinations, meaning that relaxation tends to erase discontinuities. The sharp discontinuity at cos I = 0 has already been washed away after 1 t rh, regardless of the initial amplitude of the discontinuity. We refer to appendix K.2 for a comparison of the early relaxation between a discontinuous and a smooth distribution of orbital inclinations. Note that the redistribution of prograde and retrograde orbits observed in Fig. 5 does not conflict with the conservation of the cluster's total angular momentum. Indeed, the average of cos I = Lz / L is not constrained by any global invariance, contrary to Lz . As such, even if the number of particles with positive and negative cos I (hence Lz ) changes, the conservation of the total angular momentum is ensured by a modulation of the norm of each star's angular momentum vector.  \nComparing Figs. 4 and 5 shows that orbital inclinations relax much faster than angular momenta. This is a point already raised in Rauch & Tremaine (1996) (fig. 2 therein). In particular, they showed in their section 1.4 that the long-term relaxation of E and L in spherical potential was driven by the NR theory, whereas that of the angular momentum vector L - and thus Lz - 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 signalto-noise ratio. Yet, such a projection comes at the cost of more intricate theoretical predictions, that require additional integrations along some third action. In this section, following the same approach as in T22, we first focus on in-plane relaxation, that is, relaxation occurring in ( J r , L ). Our goal is to compare quantitatively the N -body measurements with the NR prediction.  \nFig. 4. Evolution of the angular momentum DF, F ( L ), as measured in N -body simulations with N = 10 4 , with t 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.  \n<!-- image -->  \nFig. 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.  \n<!-- image -->", '4.1. N -body measurements': 'The onset of core collapse is illustrated in Fig. 6 for clusters with N = 10 5 stars. As in Fig. 3, we recover that rotation does not really impact the rate of core collapse in these linearly stable clusters. Conversely, when the clusters are unstable, such as for ( q , α ) = (1 , 0 . 5) (see Fig. 2), we find numerically (not reported in Fig. 6) that the clusters flatten and that relaxation is accelerated.  \nTo probe more precisely relaxation, let us now measure the relaxation rate, ∂ F /∂ t , in ( J r , L ). This is presented in Fig. 7 with details spelled out in Appendix E. Even in the presence of rotation, we recover the same result as in T22, namely all clusters seem to isotropise toward an in-plane DF depending only on energy. Indeed, in radially anisotropic clusters, relaxation depletes radial orbits. And the converse holds for tangentially anisotropic clusters. Importantly, the amount of rotation, α , has no significant impact on the geometry of the in-plane di ff usion. Increasing α only (very) weakly accelerates the in-plane relaxation. This is in concordance with the slightly shorter core collapse time observed in Fig. 3.  \nFig. 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 a ff ects the core contraction.  \n<!-- image -->  \nFig. 7. Relaxation rate, ∂ F /∂ t , in ( J r , L ) for various velocity anisotropies q (left to right) and rotation parameters α (top and bottom), as measured in N -body simulations. The overall amplitude of the in-plane relaxation rate depends on q , but only very weakly on α .  \n<!-- image -->', '4.2. Non-resonant prediction': 'Let us now compute the NR prediction for the relaxation rate in ( J r , L ). This prediction was numerically very intensive and required the evaluation of five embedded integrals, as demonstrated in Appendices A and H. Importantly, we stress that in the present case, the LBD (equation 5) allowed us to analytically perform the angular part of the orbit average (equation 4). This helped reducing the computation time. We refer to Appendix C.1 for details.  \nThe relaxation rate, ∂ F /∂ t in ( J r , L ), presented in Fig. 8 for various anisotropies and rotations, suggests once again an orbital reshu ffl ing toward isotropisation. In addition, increasing the rotation parameter α only very weakly increases the relaxation rate in ( J r , L ), while barely impacting the structures observed in the non-rotating case.  \nWhen comparing the N -body measurements (Fig. 7) with the NR prediction (Fig. 8), the NR theory successfully recovers the in-plane relaxation of stable rotating clusters: both figures exhibit strikingly similar structures in action space. We obtain results for the non-rotating case that are in agreement with T22. Yet, even though amplitudes in both approaches are comparable, there is still a (slight) overall prefactor mismatch. The inclusion of rotation has little impact on this mismatch. A more quantitative comparison would require performing (many) more N -body runs, as well as a more precise measurement of the relaxation rate. This would be no light undertaking, and will be the subject of future works.', '5. Out-of-plane diffusion': 'Because they rotate, the present clusters also undergo some outof-plane di ff usion, namely relative to the orbital inclinations cos I . In this section, we focus on relaxation in ( J r , cos I ), and we refer to appendix J for a similar investigation in ( L , cos I ).', '5.1. N -body measurements': 'The relaxation rates in N -body simulations are presented in Fig. 9. As expected, in the absence of rotation (top row), the  \nrelaxation is independent of cos I , and the distribution of orientations remains uniform.  \nLet us now consider the rotating clusters presented in the middle and bottom rows of Fig. 9. In radially anisotropic ( q = 1) and isotropic ( q = 0) clusters, the systems lose stars in the prograde region (cos I > 0) and gain stars in the retrograde region (cos I < 0). As one increases α (i.e. as one increases the net rotation), this trend strengthens. This is in agreement with the orbital reshu ffl ing observed in Fig. 8, where the clusters always tend to isotropise their in-plane distribution. Finally, we note that the highest di ff usion rates in inclination are observed near cos I = 0. This corresponds precisely to the location of the discontinuity of the LBD.  \nAt first glance, the tangentially anisotropic case ( q = -6), as given by the right column in Fig. 9, may seem di ff erent. Indeed, in that case a depletion of orbits for small J r is observed, whatever cos I . Conversely, the number of orbits systematically increases for large J r, whatever cos I . This is directly linked to the in-plane isotropisation of the cluster (Fig. 8). Indeed, these clusters being tangentially biased, the in-plane di ff usion occurs towards higher J r. Overall, all clusters, independently of their anisotropies, evolve towards smoother distribution of orbital inclinations.  \nAs a complement of Fig. 9, appendix J also presents the relaxation rates in ( L , cos I ). It reaches similar results as for diffusion in ( J r , cos I ). Namely, the in-plane distributions tend to isotropise and orbital inclinations di ff use so as to reach smoother distributions.  \nIn addition, we stress that due to the use of binning and finite di ff erentiation to compute ∂ F /∂ t , the expected sharp evolution occurring at cos I = 0 has been smoothed out. Furthermore, the amplitude of relaxation depends on the time interval used for finite di ff erentiation and on the bin size. While this does not drastically change the observations away from cos I = 0, any measurement in the neighborhood of cos I = 0 is tricky to perform: it should be taken with appropriate caution.  \nFig. 8. Relaxation rate, ∂ F /∂ t , in ( J 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 that depends on q , and (very) weakly on α .  \n<!-- image -->  \nFig. 9. Relaxation rate, ∂ F /∂ t , in ( J r , cos I ) for various anisotropies q (left to right) and rotation α (from top to bottom), as measured in N -body simulations. The amplitudes and structures observed depend on anisotropy, and show a reshu ffl ing towards isotropisation. Orientations redistribute toward a more a ffi ne distribution in cos I . As a result of using binning in action space and finite di ff erence in time to compute ∂ F /∂ t , the expected initially sharp relaxation occurring at cos I = 0 has been smoothed and is not visible.  \n<!-- image -->', '5.2. Non-resonant prediction': "Let us now perform the same investigation in ( J r , cos I ) using the NR theory. Obtaining a satisfying prediction near cos I = 0 required finely sampling the NR orbital integrals, as detailed in Appendix H.2. The NR prediction is illustrated in Fig. 10.  \nFirst, in the absence of rotation, the NR theory predicts a relaxation rate independent of cos I : this is a reassuring sanity  \ncheck. 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 di ff er.  \nThe di ff erence is much more striking in the presence of rotation (middle and bottom rows of Fig. 10). In both figures, the orbital reshu ffl ing still operates, and the clusters redistribute their orbits towards smoother distributions of inclinations. Yet, the  \nFig. 10. Relaxation rate, ∂ F /∂ t , in ( J r , cos I ) as in Fig. 9, but here from the NR theory (see Appendix H.2 for details). While the in-plane isotropisation and the smoothing of the inclination distribution is recovered, finer details disagree. Here, the amount of relaxation was estimated using finite di ff erentiation using separately regions with positive or negative cos I In e ff ect, this removes the δ ' (cos I ) singularity in cos I = 0, as discussed in appendix I. Importantly, the removal of this singularity does not a ff ect the prediction away from cos I = 0. In addition, because this singularity is diluted in N -body measurements - as discussed in Fig. 9 - we do not represent the NR prediction here.  \n<!-- image -->  \ndi ff usion structures predicted by the NR theory do not align with the N -body ones. Indeed, N -body simulations (Fig. 9) present 'round' structures as one goes from cos I = 0 to cos I = 1. In addition, the relaxation rate decreases as one considers orbits within to rotation plane, that is, as one considers cos I →± 1. This is in sharp contrast with the 'straighter' structures predicted by the NR theory (Fig. 10). Nevertheless, one should keep in mind that Fig. 10 cuts out the singularity at cos I = 0. The comparison to the N -body measurements near cos I = 0 is made di ffi cult by the fact that di ff usion at early time is very fast in that region, as shown by Fig. K.4. As a consequence, the measured relaxation rate may su ff er from truncation errors - due to, for example, the finite differentiation scheme - the closer one gets to cos I = 0. However, as Fig. K.4 shows, the initial discontinuity appears to have little impact on di ff usion both far away from cos I = 0 and beyond t = 0 . 1 t rh. The obtention of a more robust measure would require both a finer time di ff erence and a finer binning of action space. These two improvements would require an ensemble-average over much more realisations to reduce dispersion. A similar conclusion is reached in appendix J when considering relaxation in ( L , cos I ).  \nBefore concluding, let us elaborate further on the failure of the NR theory to model out-of-plane relaxation. For α , 0, recall that cos I = ± 1 corresponds to orbits within the rotation plane, while cos I = 0 corresponds to orbits perpendicular to the rotation plane. A test star orbiting in this rotating cluster will torque with all the other stars of the cluster. Equivalently, this means that this star will typically be subject to a torque with the total angular momentum of the rotating cluster, ⟨ L ⟩∝ α e z , the strength of which depends on the test star's orbital orientation. This is what  \nwe observe in Fig. 9. However, the NR theory only takes into account local deflections, and therefore cannot take these coherent interactions into account.  \nFollowing the work of Meiron & Kocsis (2019) (see also references therein), we anticipate for VRR to play a key role in these coherent interactions - even after factoring the truncation errors due to the finite di ff erence scheme and the singularity at cos I = 0, both discussed earlier in this section. Indeed, VRR, which is driven by persistent torques between orbital planes, might be needed to explain (i) the faster out-of-plane relaxation (section 3.3); (ii) the discrepancies between the N -body measurements and the NR predictions. Investigating quantitatively the e ffi ciency of VRR in these systems will be the topic of future investigations.", '6.1. Conclusion': "In this paper, we used the NR formalism to probe the long-term time-evolution of rotating anisotropic globular clusters. First, using N -body simulations, we showed how rotation only has a weak impact on the time of core collapse. As such, in contrast with prior research (see, e.g., Hachisu 1979), we did not observe any gravo-gyro catastrophe that could expedite core collapse (see, e.g., Einsel & Spurzem 1999). Since we introduced rotation in a di ff erent fashion, this does not contradict previous results.  \nFocusing on in-plane di ff usion, we showed how the NR prediction successfully recovers all the intricacies of the N -  \nbody measurements. This is in line with the non-rotating results from T22. Yet, although the di ff usion structures in action space closely align, there is still an overall amplitude mismatch between the N -body measurements and the kinetic prediction.  \nWe subsequently turned our interest on out-of-plane relaxation, that is, the redistribution of orbital inclinations. We pointed out the similarities and di ff erences between the N -body measurements and the NR prediction. In both approaches, we observe a systematic reshu ffl ing of orbital inclinations from overpopulated regions to underpopulated ones. As such, the distribution of inclinations get smoother through time, and this process is accelerated by rotation. On the one hand, we observed that the initial discontinuity at cos I = 0 induced a fast relaxation at small times. Surely, this strongly a ff ects any measurement of the relaxation rate in that region. On the other hand, for orbits aligned with the cluster's rotation plane - where di ff usion is much slower and measurements are much easier - the local NR theory does not match with the N -body measurements. We argue that this mismatch originates from the fact that the NR formalism is sourced by local deflections, which neglect the coherent torques between the orbital planes. Accounting for the (e ffi cient) contribution of persistent torques requires the use of VRR.", '6.2. Perspectives': 'Here, the NR formalism was used to describe the very onset of relaxation. Yet, the short duration of our simulations and our limited sample of initial conditions prevented us from fully assessing the asymptotic late time orbital distribution. Naturally, it would be interesting to push the N -body integrations further in time. While some studies carried out such long-term N -body simulations (Tiongco et al. 2020; Livernois et al. 2022), orbital inclinations were, unfortunately, not their focus. Along the same line, it would also be of interest to integrate the 3D FP equation itself. In isotropic clusters, this process has been conducted, as demonstrated for instance in the work of Vasiliev (2015). However, the scenario involving anisotropic, rotating systems has yet to be comprehensively investigated.  \nUltimately, if rotation is su ffi ciently large, it may induce a flattening of the cluster. In that case, Stäckel systems (Dejonghe & de Zeeuw 1988) could be used to model realistic flattened rotating structures. In particular, this would still ensure the integrability of the mean potential, that is the existence of explicit angle-action coordinates. Tailoring secular theory to such setups will be the topic of future works.  \nIn section 5, we pointed out how the NR theory does not reproduce the out-of-plane di ff usion observed in N -body simulations. To understand this mismatch, it would be valuable to examine the predictions of the resonant relaxation (RR) formalism (see, e.g., Hamilton et al. 2018; Fouvry et al. 2021) in the presence of rotation. While the inhomogeneous Landau equation was already investigated in non-rotating isotropic clusters (Hamilton et al. 2018; Fouvry et al. 2021), its implementation in rotating spheres with an explicit dependence on Lz should be the next step. Furthermore, to account for collective e ff ects, an explicit implementation of the Balescu-Lenard equation (Heyvaerts 2010) might prove necessary to match the details of stacked N -body measurements.  \nSimilarly, to get a better handle on the mismatch reported in section 5, one should investigate quantitatively the e ffi ciency of VRRinthese 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  \nstudied 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).  \nIt would be interesting to leverage multiple masses into the theory, so as to capture (radial or inclination) segregation effects (Meiron & Kocsis 2019). Dekel et al. (2023) argues for instance that the excess of massive galaxies found at very high redshift by the James Webb Space Telescope could originate from early feedback-free star formation within dense, possibly rotating stellar clusters made of massive stars. The fate of such clusters should be captured by the present secular theory. In particular, it could lead to the formation of intermediate mass BHs (Greene & Ho 2004; Greene et al. 2020) and possibly seeds for supermassive BHs (Kormendy & Ho 2013).', 'Data availability': 'The data underlying this article is available through reasonable request to the authors. The code, written in ju -lia (Bezanson et al. 2017), computing the NR di ff usion coe ffi -cients 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-19CE31-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 and Grant No. AST-2310362. We are grateful to M. Roule, M. Petersen, A. L. Varri and the members of the Segal collaboration for numerous suggestions during the completion of this work. We thank Stéphane Rouberol for the smooth running of the Infinity cluster, where the simulations were performed.  \nFinally, we thank our referee, Prof. Douglas Heggie, for his remarks, which contributed greatly to the improvement of this paper.", 'References': "Aarseth, S. J., Hénon, M., & Wielen, R. 1974, A&A, 37, 183 Akiyama, K. & Sugimoto, D. 1989, PASJ, 41, 991 Bar-Or, B. & Alexander, T. 2016, ApJ, 820, 129 Bellini, A. et al. 2017, ApJ, 842, 6 Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. 2017, SIAM Review, 59, 65 Bianchini, P., Varri, A. L., Bertin, G., & Zocchi, A. 2013, ApJ, 772, 67 Bianchini, P. et al. 2018, MNRAS, 481, 2125 Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton Univ. Press) Breen, P. G., Varri, A. L., & Heggie, D. C. 2017, MNRAS, 471, 2778 Casertano, S. & Hut, P. 1985, ApJ, 298, 80 Chandrasekhar, S. 1943, ApJ, 97, 255 Cohn, H. 1979, ApJ, 234, 1036 Dejonghe, H. 1986, Phys. Rep., 133, 217 Dejonghe, H. 1987, MNRAS, 224, 13 Dejonghe, H. & de Zeeuw, T. 1988, ApJ, 333, 90 Dekel, A., Sarkar, K. C., Birnboim, Y., Mandelker, N., & Li, Z. 2023, MNRAS, 523, 3201 Eilon, E., Kupi, G., & Alexander, T. 2009, ApJ, 698, 641 Einsel, C. & Spurzem, R. 1999, MNRAS, 302, 81 Ernst, A., Glaschke, P., Fiestas, J., Just, A., & Spurzem, R. 2007, MNRAS, 377,  \n- 465  \nFabricius, M. H., Noyola, E., Rukdee, S., et al. 2014, ApJL, 787, L26  \nFerraro, F. R., Mucciarelli, A., Lanzoni, B., et al. 2018, ApJ, 860, 50 Fiestas, J. & Spurzem, R. 2010, MNRAS, 405, 194 Foote, H. R., Generozov, A., & Madigan, A.-M. 2020, ApJ, 890, 175 Fouvry, J.-B., Bar-Or, B., & Chavanis, P.-H. 2019, ApJ, 883, 161 Fouvry, J.-B., Hamilton, C., Rozier, S., & Pichon, C. 2021, MNRAS, 508, 2210 Geyer, E. H., Hopp, U., & Nelles, B. 1983, A&A, 125, 359 Giersz, M. & Heggie, D. C. 1994, MNRAS, 268, 257 Ginat, Y. B., Panamarev, T., Kocsis, B., & Perets, H. B. 2023, MNRAS, 525,  \n4202 Greene, J. E. & Ho, L. C. 2004, ApJ, 610, 722 Greene, J. E., Strader, J., & Ho, L. C. 2020, ARA&A, 58, 257 Gruzinov, A., Levin, Y., & Zhu, J. 2020, ApJ, 905, 11 Hachisu, I. 1979, PASJ, 31, 523 Hachisu, I. 1982, PASJ, 34, 313 Hamilton, C., Fouvry, J.-B., Binney, J., & Pichon, C. 2018, MNRAS, 481, 2041 Heggie, D. & Hut, P. 2003, The Gravitational Million-Body Problem (Cambridge University Press) 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, J. 2010, MNRAS, 407, 355 Inagaki, S. & Hachisu, I. 1978, PASJ, 30, 39 Jindal, A., Webb, J. J., & Bovy, J. 2019, MNRAS, 487, 3693 Kamann, S., Husser, T.-O., Dreizler, S., et al. 2018, MNRAS, 473, 5591 Kamlah, A. W. H., Spurzem, R., Berczik, P., et al. 2022, MNRAS, 516, 3266 Kim, E., Einsel, C., Lee, H. M., Spurzem, R., & Lee, M. G. 2002, MNRAS, 334, 310 Kim, E., Lee, H. M., & Spurzem, R. 2004, MNRAS, 351, 220 King, I. R. 1966, AJ, 71, 64 Kocsis, B. & Tremaine, S. 2011, MNRAS, 412, 187 Kocsis, B. & Tremaine, S. 2015, MNRAS, 448, 3265 Kontizas, E., Kontizas, M., Sedmak, G., & Smareglia, R. 1989, AJ, 98, 590 Kormendy, J. & Ho, L. C. 2013, ARA&A, 51, 511 Lagoute, C. & Longaretti, P. Y. 1996, A&A, 308, 441 Lanzoni, B., Ferraro, F. R., Mucciarelli, A., et al. 2018, ApJ, 865, 11 Livernois, A. R., Vesperini, E., Varri, A. L., Hong, J., & Tiongco, M. 2022, MNRAS, 512, 2584 Lynden-Bell, D. 1960, MNRAS, 120, 204 Lynden-Bell, D. & Wood, R. 1968, MNRAS, 138, 495 Magnan, N., Fouvry, J.-B., Pichon, C., & Chavanis, P.-H. 2022, MNRAS, 514, 3452 McLaughlin, D. E. & van der Marel, R. P. 2005, ApJs, 161, 304 Meiron, Y. & Kocsis, B. 2019, ApJ, 878, 138 Miocchi, P., Lanzoni, B., Ferraro, F. R., et al. 2013, ApJ, 774, 151 Rauch, K. P. & Tremaine, S. 1996, New. Astron., 1, 149 Risken, H. 1996, The Fokker-Planck Equation, 2nd edn. (Berlin: Springer), gst: Rosenbluth, M. N., MacDonald, W. M., & Judd, D. L. 1957, Phys. Rev., 107, 1 Rozier, S., Fouvry, J.-B., Breen, P. G., et al. 2019, MNRAS, 487, 711 Sollima, A., Baumgardt, H., & Hilker, M. 2019, MNRAS, 485, 1460 Spitzer, L. 1975, in Dynamics of Stellar Systems, Vol. 69, 3 Stetson, P. B., Pancino, E., Zocchi, A., Sanna, N., & Monelli, M. 2019, MNRAS, 485, 3042 Szölgyén, Á. & Kocsis, B. 2018, PRL, 121, 101101 Szölgyén, Á., Máthé, G., & Kocsis, B. 2021, ApJ, 919, 140 Szölgyén, Á., Meiron, Y., & Kocsis, B. 2019, ApJ, 887, 123 Tep, K., Fouvry, J.-B., & Pichon, C. 2022, MNRAS, 514, 875 Tiongco, M., Collier, A., & Varri, A. L. 2021, MNRAS, 506, 4488 Tiongco, M., Vesperini, E., & Varri, A. L. 2016, MNRAS, 461, 402 Tiongco, M., Vesperini, E., & Varri, A. L. 2020, in Star Clusters: From the Milky Way to the Early Universe, Vol. 351, 524-527 Tiongco, M., Vesperini, E., & Varri, A. L. 2022, MNRAS, 512, 1584 Trager, S. C., King, I. R., & Djorgovski, S. 1995, AJ, 109, 218 Varri, A. L. & Bertin, G. 2012, A&A, 540, A94 Vasiliev, E. 2015, MNRAS, 446, 3150 Wang, L., Spurzem, R., Aarseth, S., et al. 2015, MNRAS, 450, 4070 Watkins, L. L., van der Marel, R. P., Bellini, A., & Anderson, J. 2015, ApJ, 803, 29  \nWhite, R. E. & Shawl, S. J. 1987, ApJ, 317, 246  \nWilson, C. P. 1975, AJ, 80, 175  \nFig. A.2. Illustration of the frames used to compute the argument of the background DF, F rot( E ' , L ' , L ' z ). 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.  \n<!-- image -->  \nFig. A.1. Tailored frame used to compute the parallel and perpendicular local velocity deflections in equations (A.1). By construction, the test star's angular momentum L is along the axis e 2.  \n<!-- image -->", '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.', 'Appendix A.1: Local velocity coefficients': "The starting point are equations (2). We consider this relation within the frame from Fig. A.1. The local velocity deflections are then given by  \n⟨ ∆ v ∥ ⟩ = ⟨ ∆ v 1 ⟩ , (A.1a) ⟨ ( ∆ v ∥ ) 2 ⟩ = ⟨ ( ∆ v 1) 2 ⟩ , (A.1b)  \n⟨ ( ∆ v ⊥ ) 2 ⟩ = ⟨ ( ∆ v 2) 2 ⟩ + ⟨ ( ∆ v 3) 2 ⟩ , (A.1c)  \nwhere the coordinate 1 is the z coordinate, parallel to v . The coordinates 2 and 3 are defined such that the projection of r on the ( Oxy ) plane is along e 3 (Fig. A.1). We now introduce w = v -v ' and change the integration variables from v ' to w in equation (2). We get  \n⟨ ∆ vi ⟩ = -2 A Z d w wi w 3 F rot , (A.2a)  \n⟨ ∆ vi ∆ vj ⟩ = A Z d w w 2 δ i j -wiwj w 3 F rot , (A.2b)  \nwith the shortened notation F rot = F rot( r , v ' = v -w ). Using these equations allows us to compute the needed coe ffi cients from equation (A.1) through  \n               ⟨ ∆ v 1 ⟩ ⟨ ( ∆ v 1) 2 ⟩ ⟨ ( ∆ v 2) 2 ⟩ ⟨ ( ∆ v 3) 2 ⟩                = A Z d w w 3                   -2 w 1 w 2 -w 2 1 w 2 -w 2 2 w 2 -w 2 3                   F rot . (A.3)  \nWe now define the spherical coordinates associated with this frame (Fig. A.1)  \nw 1 = w cos ϑ ; w 2 = w sin ϑ cos ϕ ; w 3 = w sin ϑ sin ϕ. (A.4)  \nInjecting equations (A.4) into equations (A.3), we obtain  \n               ⟨ ∆ v 1 ⟩ ⟨ ( ∆ v 1) 2 ⟩ ⟨ ( ∆ v 2) 2 ⟩ ⟨ ( ∆ v 3) 2 ⟩                = A Z d w d ϑ d ϕ sin ϑ                 -2 cos ϑ w sin 2 ϑ w (1 -sin 2 ϑ cos 2 ϕ ) w (1 -sin 2 ϑ sin 2 ϕ )                 F rot . (A.5)  \nHere, the background distribution, F rot, is to be evaluated in F rot( r , v ' ) = F rot( E ' , L ' , L ' z ). We compute these arguments in appendix A.2. Finally injecting equation (A.5) into equation (A.1), we obtain our main result, namely equation (3), as given in the main text.", 'Appendix A.2: Arguments of the background distribution': "Let us now compute the arguments E ' , L ' and L ' z for the background distribution function, F rot( E ' , L ' , L ' z ), in equation (2).  \nWe start with E ' . It reads  \nE ' = ψ ( r ) + v ' 2 2 = ψ ( r ) + ( v -w ) 2 2 = E + w 2 2 -wv cos ϑ, (A.6)  \nsince v = v e 1 (see Fig. A.1). This form allows us to get an upper bound on the w -integral of equation (A.5), above which E ' > 0. As soon as E ' > 0, the background star is unbound and its DF vanishes.  \nThe computation of L ' and L ' z is best achieved by using a di ff erent coordinate system, as illustrated in Figs. A.1 and A.2. More precisely, the frame (1 , 2 , 3) is related to the ( X , Y , Z ) one through the relations  \n          e 1 e 2 e 3           =           v r v cos φ -v t v sin φ v r v sin φ + v t v cos φ 0 0 0 1 v r v sin φ + v t v cos φ v t v sin φ -v r v cos φ 0                     e X e Y e Z           , (A.7)  \nwith v r (resp. v t) the radial (resp. tangential) component of the test star's velocity. We can finally relate the frame ( X , Y , Z ) to the frame ( x , y , z ) through the relations  \n          e X e Y e Z           =           cos I 0 sin I 0 1 0 -sin I 0 cos I                     e x e y e z           . (A.8)  \nWe can now compute L , the norm of L . As it is a norm, we can choose any frame to compute it. The frame (1 , 2 , 3) is a convenient choice as this frame is used to compute the integrands of equations (A.3). We write  \nL ' 2 = ( r 1 v ' 2 -r 2 v ' 1 ) 2 + ( r 1 v ' 3 -r 3 v ' 1 ) 2 + ( r 2 v ' 3 -r 3 v ' 2 ) 2 , (A.9)  \nwhere ( r 1 , r 2 , r 3) are the coordinates of r in the frame (1 , 2 , 3). Using equations (A.7) and (A.8), we obtain  \nr 1 = v r v r ; r 2 = 0 ; r 3 = v t v r . (A.10)  \nArticle number, page 11 of 18  \nTherefore, one gets", 'Appendix B: 3D diffusion coefficients': "L ' = r s v ' 2 2 + v r v v ' 3 -v t v v ' 1 ! 2 (A.11) = r s ( w sin ϑ cos ϕ ) 2 + v t + v r v w sin ϑ sin ϕ -v t v w cos ϑ ! 2 ,  \ngiven that v ' = v -w .  \nLet us finally compute L ' z . In the ( x , y , z ) frame, it reads  \nL ' z = xv ' y -yv ' x . (A.12)  \nFrom equation (A.8), the positions read  \nx = X cos I = r cos φ cos I , (A.13a)  \ny = Y = r sin φ, (A.13b)  \nz = X sin I = r cos φ sin I , (A.13c)  \nusing X = r cos φ and Y = r sin φ . We then obtain the velocities  \nvx = v r cos φ cos I -v t sin φ cos I , (A.14a)  \nvy = v r sin φ + v t cos φ, (A.14b)  \nvz = v r cos φ sin I -v t sin φ sin I , (A.14c)  \nby taking the time derivative of equations (A.13). Let us now compute the background velocity components. First, using the left panel of Fig. A.2, we have  \n          v ' x v ' y v ' z           =           cos I 0 -sin I 0 1 0 sin I 0 cos I                     v ' X v ' Y v ' Z           . (A.15)  \nSecond, using the relation in equations (A.7) (Fig. A.1 and right panel of Fig. A.2), we have  \n          v ' X v ' Y v ' Z           =           v r v cos φ -v t v sin φ 0 v r v sin φ + v t v cos φ v r v sin φ + v t v cos φ 0 v t v sin φ -v r v cos φ 0 1 0                     v ' 1 v ' 2 v ' 3           . (A.16)  \nOverall, we can rewrite equation (A.12) as  \nL ' z = r v ' 1 v t v -v ' 3 v r v ! cos I + rv ' 2 sin φ sin I = r v t -v t v w cos ϑ + v r v w sin ϑ sin ϕ ! cos I -rw sin ϑ cos ϕ sin φ sin I . (A.17)  \nFinally, using the Cauchy-Schwarz inequality, one can check that  \n| L ' z | ≤ L ' q cos 2 I + sin 2 φ sin 2 I , (A.18)  \nwhich ensures that | L ' z | ≤ L ' .  \nEquations (A.6), (A.11) and (A.17) are the main results of this section. They provide us with explicit expressions for the arguments at which to evaluate F rot = F rot( E ' , L ' , L ' z ) in equation (3).  \nThe local in-plane di ff usion coe ffi cients, that is, the di ff usion coe ffi cients in E and L , are given in appendix C of Bar-Or & Alexander (2016). We reproduce them here for completeness. They read  \n⟨ ∆ E ⟩ = 1 2 ⟨ ( ∆ v ∥ ) 2 ⟩ + 1 2 ⟨ ( ∆ v ⊥ ) 2 ⟩ + v ⟨ ∆ v ∥ ⟩ , (B.1a)  \n⟨ ( ∆ E ) 2 ⟩ = v 2 ⟨ ( ∆ v ∥ ) 2 ⟩ , (B.1b)  \n⟨ ∆ L ⟩ = r v t v ⟨ ∆ v ∥ ⟩ + r 2 4 L ⟨ ( ∆ v ⊥ ) 2 ⟩ , (B.1c)  \n⟨ ( ∆ L 2 ) ⟩ = r 2 v 2 t v 2 ⟨ ( ∆ v ∥ ) 2 ⟩ + r 2 2 v 2 r v 2 ⟨ ( ∆ v ⊥ ) 2 ⟩ , (B.1d)  \n⟨ ∆ E ∆ L ⟩ = L ⟨ ( ∆ v ∥ ) 2 ⟩ . (B.1e)  \nIn order to compute the 3D FP equation (1), we also need the local di ff usion coe ffi cients involving Lz = L · e z . Performing a firstorder variation of L = r × v , we can write  \n∆ L = ∆ v ∥ v -v r v t ∆ v 3 v ! L + L ∆ v 2 v v r v t e 3 -ˆ v ! , (B.2a)  \ne 3 = v × L | v × L | = v ˆ r -v r ˆ v v t , (B.2b)  \nwhere ( e 1 , e 2 , e 3) is the frame illustrated in Fig. A.1. We can then write  \ne 3 · e z = sin I v r v sin φ + v t v cos φ ! . (B.3)  \nAs a result, it follows that  \nv r v e 3 · e z -vzv t v 2 = sin φ sin I , (B.4)  \nwith φ given in Fig. (A.2) and vz by equation (A.14c). Overall, we finally have  \n⟨ ∆ Lz ⟩ = Lz v ⟨ ∆ v ∥ ⟩ , (B.5a)  \n⟨ ( ∆ L 2 z ⟩ = Lz L ! 2 L 2 v 2 ⟨ ( ∆ v ∥ ) 2 ⟩ + 1 2 r 2 v 2 r v 2 ⟨ ( ∆ v ⊥ ) 2 ⟩ ! (B.5b) + r 2 sin 2 φ 2 1 -L 2 z L 2 ! ⟨ ( ∆ v ⊥ ) 2 ⟩ ,  \n⟨ ∆ E ∆ Lz ⟩ = Lz ⟨ ( ∆ v ∥ ) 2 ⟩ , (B.5c)  \n⟨ ∆ L ∆ Lz ⟩ = Lz L L 2 v 2 ⟨ ( ∆ v ∥ ) 2 ⟩ + 1 2 r 2 v 2 r v 2 ⟨ ( ∆ v ⊥ ) 2 ⟩ ! . (B.5d)  \nEquations (B.1) and (B.5) are the main results of this section. They allow us to compute the local di ff usion coe ffi cients in ( E , L , Lz ) from the local velocity di ff usion coe ffi cients given by equation (3).", 'Appendix C: Orbit average': 'The computation of the global di ff usion coe ffi cients in the FP equation (1) requires that we orbit average the local di ff usion coe ffi cients 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.  \nTable C.1. Explicit values of A 1 and A 2 (equation C.2a) for the LBD DF (equation 5).', 'Appendix C.1: Angular orbit average': "We start from the generic formula of equation (4). There, the average with respect to φ must be performed with care since L ' z depends on φ (equation A.17). Fortunately, taking advantage of the particular structure of the di ff usion coe ffi cients in E , L and Lz (see equations B.1 and B.5), we only need to compute the following four integrals over φ  \nZ d φ 2 π ⟨ ∆ v ∥ ⟩ , Z d φ 2 π ⟨ ( ∆ v ∥ ) 2 ⟩ , Z d φ 2 π ⟨ ( ∆ v ⊥ ) 2 ⟩ , Z d φ 2 π cos 2 φ ⟨ ( ∆ v ⊥ ) 2 ⟩ , (C.1)  \nwhere we used sin 2 φ = 1 2 (1 -cos 2 φ ) for convenience.  \nIn equation (3), the only dependence with respect to φ is in the dependence with respect to L ' z in F rot. In the particular case of the LBD (equation 5), we can perform explicitly the average over φ . Computing the integrals (C.1) only requires the evaluation of the two non-trivial integrals  \nA 1 = Z 2 π d φ 2 π sgn( L ' z / L ' ) ,  \nA 2 = Z 2 π 0 d φ 2 π cos 2 φ sgn( L ' z / L ' ) ,  \n0 (C.2a) (C.2b)  \nwith L ' z following from equation (A.17). Introducing ζ = -( v ' 1 v t / v -v ' 3 v r / v ) cos I and µ = ζ/ ( v ' 2 sin I ), these two integrals can be explicitly computed. Their expressions are gathered in Table C.1. This is the main result of this section.  \nFrom 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.", 'Appendix C.2: Radial orbit average': 'Having performed the angle average in Appendix C.1, equation (4) takes the form  \nDX = Ω r π Z r a r p d r | v r | ⟨ ∆ X ⟩ φ ( r ) , (C.3)  \nwhere ⟨·⟩ φ stands for the angle average. As described in appendix F2 of T22, equation (C.3) can be e ffi ciently performed by changing variable from the radius to an e ff ective anomaly. In particular, this removes the integrable singularity 1 / | v r | at pericentre and apocentre. We follow the exact same approach here, and perform the radial average using 50 sampling nodes.', 'Appendix D: Fokker-Planck and coordinates': 'Depending on the quantities we wish to investigate, we may wish to change of coordinates system. Fortunately, it is possible to transform one FP description into another under a change of variable. We describe a few relevant examples in this appendix.', 'Appendix D.1: Generic change of coordinates': 'Given some coordinates x , the FP equation generically reads (Risken 1996)  \n∂ F ( x , t ) ∂ t = -∂ ∂ x · F x ( x ) (D.1) = -∂ ∂ x · " D x ( x ) F ( x ) -1 2 ∂ ∂ x · D xx ( x ) F ( x ) !# ,  \nwith the di ff usion coe ffi cients  \nD x ( x ) =           Dx 1 Dx 2 Dx 3           ; D xx ( x ) =           Dx 1 x 1 Dx 1 x 2 Dx 1 x 3 Dx 1 x 2 Dx 2 x 2 Dx 2 x 3 Dx 1 x 3 Dx 2 x 3 Dx 3 x 3           , (D.2)  \nand the DF in x space, F ( x ). In our case, depending on the context, x = ( x 1 , x 2 , x 3) may either stand for ( E , L , Lz ), ( J r , L , Lz ) or ( J r , L , cos I ).  \nFollowing Risken (1996) and Bar-Or & Alexander (2016), one can easily rewrite the FP equation within some new coordinates x \' ( x ). The new di ff usion coe ffi cients read  \nDx \' l = X k ∂ x \' l ∂ xk Dx k + 1 2 X k , r ∂ 2 x \' l ∂ xr ∂ xk Dx r xk , (D.3a)  \nDx \' l x \' m = X k , r ∂ x \' l ∂ xr ∂ x \' m ∂ xk Dx r xk . (D.3b)  \nThese coe ffi cients source a FP equation in x \' -space for the DF, F \' ( x \' ), reading  \nF \' ( x \' ) = GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> ∂ x ∂ x \' GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> F ( x ) , (D.4)  \nwith | ∂ x /∂ x \' | the inverse Jacobian of the coordinate transform.  \nAppendix D.2: From ( E , L , Lz ) to ( J r , L , Lz )  \nFollowing section 2, we have at our disposal di ff usion coe ffi cients within the coordinates ( E , L , Lz ). Owing to equations (D.3), the di ff usion coe ffi cients in the action space ( J r[ E , L ] , L , Lz ) are easily computed. This is already detailed in appendix F3 of T22. Ultimately, we have at our disposal the diffusion coe ffi cients  \nD 1( J ) =           DJ r DL DL z           ; D 2( J ) =           DJ r J r DJ r L DJ r Lz DJ r L DLL DLLz DJ r Lz DLLz DLzLz           , (D.5)  \nwhich source the FP evolution of F rot( J r , L , Lz ), the cluster\'s DF in ( J r , L , Lz ).  \nAppendix D.3: From ( J r , L , Lz ) to ( J r , L , cos I )  \nLet us now change of coordinates from J = ( J r , L , Lz ) to J c = ( J r , L , cos I = Lz / L ). Applying equations (D.3), we obtain the new di ff usion coe ffi cients  \nD cos I = -cos I DL L + DL z L + cos I DLL L 2 -DLL z L 2 , (D.6a)  \nDJ r cos I = -cos I DJ r L L + DJ r Lz L = 0 , (D.6b)  \nDL cos I = -cos I DLL L + DLL z L = 0 , (D.6c)  \nD cos I cos I = cos 2 I DLL L 2 -2 cos I DLL z L 2 + DL z Lz L 2 , (D.6d)  \nTable E.1. 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 r , L ) domain in NJ r × NL uniform bins within the region J min r ≤ J r ≤ J max r (similarly for L ). We use a similar approach for ( J r , cos I ) and ( L , cos I ). All quantities are in physical units G = M = b = 1, if not stated otherwise.  \nwhile the other coe ffi cients stay unchanged. Importantly, we point out that DJ r cos I = DL cos I = 0. This comes from the relations DEL z = cos I DEL and DLL z = cos I DLL , which are inferred from equations (B.1) and (B.5). Hence, we are left with the di ff usion coe ffi cients  \nD 1( J c) =           DJ r DL D cos I           ; D 2( J c) =           DJ r J r DJ r L 0 DJ r L DLL 0 0 0 D cos I cos I           . (D.7)  \nOverall, they drive the FP evolution of F ( J c) = LF rot( J ), the cluster\'s DF in ( J r , L , cos I ).', 'Appendix E: N -body simulations': "The simulations presented throughout the main text were performed using the direct N -body code NBODY6 ++ GPU (Wang et al. 2015), version 4.1. We used the exact same run parameters as in appendix G of T22, and drew the initial conditions using P lummer P lus . py (Breen et al. 2017). We also refer to the aforementioned appendix for the N -body measurements, and to Table E.1 for all our binning parameters.  \nEach N -body realisation was composed of N = 10 5 (resp. N = 10 4 ) stars and integrated up to t max = 1000 HU (resp. t max = 4000 HU) with a dump every 1HU. On a 40-core CPU node with a single V100 GPU, one simulation typically required ∼ 24 h (resp. ∼ 12 h) of computation. Ensemble averages were performed over 50 independent runs.  \nIn practice, the main di ffi culties in the N -body measurements are (i) the estimation of the instantaneous potential - necessary to compute the instantaneous radial action, J r; (ii) the determination of the bins' size in action space; (iii) the estimation of the relaxation rate ∂ F /∂ t via finite di ff erences. This is especially true for highly tangentially anisotropic clusters (e.g., q = -6), where stars are closely stacked near the J r = 0 axis. For these three difficulties, we use the same approach as in T22.", 'Appendix F: Sphericity': "In order to track the clusters' sphericity in Fig. 2, we introduce the 3D inertia-like matrix  \nI = N X k = 1 ρ 2 k r T k r k I -r k r T k !, N X k = 1 ρ 2 k , (F.1)  \nwith r k the location of the k -th particle, ρ k its local density (see Casertano & Hut 1985) and I the 3D identity matrix. In that definition, the extra ρ 2 k factors enhance the contributions from the regions close to the centre.  \nThe matrix I is obviously symmetric. It is also semi-definite positive since for any y ∈ R 3 , we have  \ny T I y = X k ρ 2 k r T k r k y T y -y T r k r T k y !, X k ρ 2 k = X k ρ 2 k | r k | 2 | y 2 | 2 -| r k · y | 2 !, X k ρ 2 k ≥ 0 , (F.2)  \nfollowing Cauchy-Schwarz's inequality. As a result, I has three positive eigenvalues, { λ i } i , which encapsulate the cluster's sphericity. Indeed, spherically symmetric clusters have all their eigenvalues equal.  \nWe generically define the cluster's sphericity via h = λ min /λ max, which we estimated from N -body simulations. To reduce shot noise in that measurement, we averaged h over realisations as follows. First, we computed the elementary symmetric polynomials α = λ 1 + λ 2 + λ 3, β = λ 1 λ 2 + λ 1 λ 3 + λ 2 λ 3 and γ = λ 1 λ 2 λ 3 for every cluster and every timestep. We then averaged the values of ( α, β, γ ) over all realisations. From these values, we estimated the eigenvalues λ i as the three (positive) roots of the polynomial λ 3 -⟨ α ⟩ λ 2 + ⟨ β ⟩ λ -⟨ γ ⟩ . 4  \nIn Fig. 2, we illustrate the evolution of the sphericity, h , for various amounts of anisotropy and rotation. We typically find h ≃ 0 . 996, that is, clusters remain reasonably spherically symmetric throughout their evolution. As as visual check of this conclusion, we represent in Fig. F.1 the late time stellar distribution of a rotating isotropic globular cluster with N = 10 4 . Even close to core collapse, the cluster remains spherically symmetric.", '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  \nF rot( E , Lz ) ∝ (e -β E -1) e -β Ω 0 Lz , (G.1)  \nwhere F rot = F rot( r , v ) = F rot( E , Lz ) is the DF in ( r , v ). In this DF, β = 1 /σ 2 c is an inverse temperature with σ c the central velocity dispersion, while ω 0 = √ 9 / (4 π Gn c) Ω 0 is a rotation parameter with n c the central density and Ω 0 the angular velocity in the cluster's centre (Lagoute & Longaretti 1996). In particular, β is related to the so-called King parameter (King 1966), defined by W 0 = -β ( ψ -ψ t) where ψ t is the potential at the cluster's edge. The rotating King models from Einsel & Spurzem (1999) can therefore be parametrised by ( W 0 , ω 0). In practice, Einsel & Spurzem (1999) integrated these models up to core collapse using a FP scheme in ( E , Lz ) for di ff erent rotation parameters, ω 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 ω 0, the cluster's mean potential also varies.  \nFig. F.1. Snapshots of the distribution of stars in a rotating ( α = 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 HU. As time evolves, the cluster contracts and the central density increases. Notwithstanding, the cluster retains its initial spherical symmetry.", 'Appendix H: 2D Fokker-Planck equations': 'In practice, it is convenient to study long-term relaxation through two-dimensional projections of action space. To do so, we must integrate over one coordinate the 3D FP equation expressed either in ( J r , L , Lz ) (appendix D.2) or in ( J r , L , cos I ) (appendix D.3). In this appendix, we compute the 2D FP equations in ( J r , L ), ( J r , cos I ) and ( L , cos I ) that respectively drive the evolution of the DFs  \nF ( J r , L ) = Z L -L d Lz F rot( J r , L , Lz ) , (H.1a)  \nF ( J r , cos I ) = Z + ∞ 0 d L F ( J r , L , cos I ) , (H.1b)  \nF ( L , cos I ) = Z + ∞ 0 d J r F ( J r , L , cos I ) . (H.1c)  \nAppendix H.1: Equation in ( J r , L )  \nWe start from equation (1), which describes di ff usion in ( J r , L , Lz ), and write the flux divergence as  \n∂ F rot ∂ t = -∂ F J r ∂ J r + ∂ F L ∂ L + ∂ F Lz ∂ Lz ! . (H.2)  \nIntegrating over Lz yields  \nZ L -L d Lz ∂ F J r ∂ J r = ∂ ∂ J r Z L -L d Lz F J r , (H.3a)  \nZ L -L d Lz ∂ F L ∂ L = -F L ( Lz = L ) -F L ( Lz = -L ) + ∂ ∂ L Z L -L d Lz F L , (H.3b)  \nZ L -L d Lz ∂ F Lz ∂ Lz = F Lz ( Lz = L ) - F Lz ( Lz = -L ) . (H.3c)  \nIn addition, because the flux cannot exit action space, we also have F Lz ( Lz = ± L ) = ±F L ( Lz = ± L ). As a result, all the nonintegral terms in equations (H.3) cancel one another. We are then left with the 2D equation  \n∂ F ( J r , L , t ) ∂ t = -∂ ∂ ( J r , L ) · Z L -L d Lz        F J r ( J r , L , Lz ) F L ( J r , L , Lz )        . (H.4)  \nIt is this rewriting that is used in section 4.2. To obtain Fig. 8, the Lz -integrals in equations (H.3) are sampled with 50 nodes. The ( w , ϑ, ϕ )-integrals of equations (3) are sampled with 100 w -nodes, 100 ϑ -nodes and 200 ϕ -nodes.', 'Appendix H.2: Equation in ( J r , cos I )': 'We start from the di ff usion coe ffi cients in equation (D.7), which describe di ff usion in ( J r , L , cos I ), for the DF, F = LF rot. We write the flux divergence as  \n∂ F ∂ t = -∂ F J r ∂ J r + ∂ F L ∂ L + ∂ F cos I ∂ cos I ! . (H.5)  \nIntegrating over L yields  \nZ ∞ 0 d L ∂ F J r ∂ J r = ∂ ∂ J r Z + ∞ 0 d L F J r , (H.6a)  \nZ ∞ 0 d L ∂ F L ∂ L = F L ( L =+ ∞ ) - F L ( L = 0) , (H.6b)  \nZ ∞ 0 d L ∂ F cos I ∂ cos I = ∂ ∂ cos I Z + ∞ 0 d L F cos I . (H.6c)  \nBecause the flux cannot leave action space, equation (H.6b) vanishes. Then, equation (H.5) readily reduces to a 2D FP equation in ( J r , cos I ), which we use in section 5.2. To obtain Fig. 10, the L -integrals in equations (H.6) are sampled with 50 nodes for 0 ≤ L ≤ 3 L 0, where L 0 = √ GMb is the typical action. The ( w , ϑ, ϕ )-integrals of equations (3) are sampled with 100 w -nodes, 100 ϑ -nodes and 800 ϕ -nodes.', 'Appendix H.3: Equation in ( L , cos I )': "As in appendix H.2, we start from the di ff usion coe ffi cients in equation (D.7). Integrating equation (H.5) over J r yields  \nZ ∞ 0 d J r ∂ F J r ∂ J r = F J r ( J r =+ ∞ ) - F J r ( J r = 0) , (H.7a)  \nZ ∞ 0 d J r ∂ F L ∂ L = ∂ ∂ L Z + ∞ 0 d J r F L , (H.7b)  \nZ ∞ 0 d J r ∂ F cos I ∂ cos I = ∂ ∂ cos I Z + ∞ 0 d J r F cos I . (H.7c)  \nOnce again, because the flux cannot leave action space, equation (H.7a) vanishes. Then, equation (H.5) reduces to a 2D FP equation in ( L , cos I ), which we use in appendix J. To obtain Fig. J.2, the J r-integrals in equations (H.7) are sampled with 50 nodes for 0 ≤ J r ≤ 10 L 0. The ( w , ϑ, ϕ )-integrals of equations (3) are sampled with 100 w -nodes, 100 ϑ -nodes and 1600 ϕ -nodes.  \nFig. I.1. Representation for a Plummer cluster ( q , α ) = (0 , 0 . 25) of the cos I terms, D cos I F rot (on the top left) and D cos I cos I F rot (on the top right), respectively with the first (on the bottom left) and second (on the bottom right) derivative with respect to cos I . In particular, the bottom right panel displays an approximate δ ' (cos I ) behaviour - the amplitude of the jumps tends to infinity as one reduces the step in the finite di ff erentiation scheme - which is the expected theoretical relaxation induced by the sgn discontinuity of F rot at t = 0.  \n<!-- image -->", 'Appendix I: Evaluating the discontinuity at cos I = 0': "In this appendix, we highlight the sharp discontinuity of the relaxation rate that occurs at cos I = 0 as a result of the discontinuous LBD distribution of orientations (equation 5). This singularity has been removed from Figures 10 and J.2 (see appendices H.2 and H.3). Because equation (5) is discontinuous, we expect that taking the first and second derivatives with respect to cos I will yield a theoretical δ ' (cos I ) behaviour around cos I = 0. To that aim, we show in Figure I.1 the quantities D cos I F rot and D cos I cos I F rot, respectively with their first and second derivatives with respect to cos I . In particular, ∂ 2 ( D cos I cos I F rot) /∂ cos I 2 has a δ ' (cos I ) component near cos I = 0. In practice, this is smoothed by the finite di ff erentiation used here. Indeed, reducing the finite di ff erentiation step sharpens the discontinuity to higher and higher values, hence converging to the true δ ' (cos I ) behaviour.", 'Appendix J: 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. J.1, we illustrate the N -body measurements while Fig. J.2 presents the associated NR predictions.  \nFrom these two figures, we can globally make the same observations as in section 5, which considered di ff usion in ( J r , cos I ). Indeed, in rotating clusters, di ff usion reshu ffl es orbits towards smoother distributions of inclinations. Nevertheless, the NRprediction 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 di ff ers from the N -body simulations, where the relaxation rate decreases away from cos I = 0.', 'Appendix K: 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  \nF rot( J r , L , Lz ) = F tot( J r , L )(1 + α g [ Lz / L ]) , (K.1)  \nwhere g [cos I ] is an odd function with g (1) = 1. The LBD corresponds to g = sgn. To approximate smoothly the LBD, we consider the sequence of functions  \nga ( x ) = erf( ax ) / erf( a ) . (K.2)  \nAs illustrated in Fig. K.1, this ensures that g 0( x ) = x and g ∞ ( x ) = sgn( x ).', 'Appendix K.1: The cos I coordinate': 'To probe the possible presence of discontinuities and singularities, let us first compute the cos I component of the 3D flux in ( J r , L , cos I ), as introduced in appendix D.3. The dependence of this flux with respect to a is illustrated in Fig. K.2.  \nIn this figure, for any smooth ga , we observe a 1 / L divergence of the flux as L → 0. In a nutshell, when working with cos I , we su ff er from a coordinate singularity, and the NR prediction cannot be applied to a DF with a smooth rotation function. Yet, when a →∞ , the flux converges towards the LBD flux pointwise, and this flux does not diverge for L → 0. Phrased differently, in the particular case of the LBD sign function, we can make a meaningful and well-posed NR prediction for the di ff usion in cos I . In practice, this vanishing of the divergence stems from the cancellation of the derivative of 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 (D.6) do not combine into an integrable quantity.  \nTo further stress that this divergence originate from coordinate singularities, let us now produce NR predictions by describing relaxation in ( J r , L , Lz ) rather than in ( J r , L , cos I ). To do so, we define the DF in Lz as  \nF ( Lz ) = Z + ∞ 0 d J r Z + ∞ | Lz | d L F rot( J r , L , Lz ) . (K.3)  \nIntegrating the 3D FP equation (1) over J r and L yields  \n∂ F ( Lz ) ∂ t = -∂ ∂ Lz GLYPH<16> DL z F [ Lz ] GLYPH<17> + 1 2 ∂ 2 ∂ L 2 z GLYPH<16> DL z Lz F [ Lz ] GLYPH<17> . (K.4)  \nHere, the 1D di ff usion coe ffi cients in Lz are given by  \nDL z = Z + ∞ 0 d J r Z + ∞ | Lz | d L DL z F rot( J r , L | Lz ) , (K.5a)  \nDL z Lz = Z + ∞ 0 d J r Z + ∞ | Lz | d L DL z Lz F rot( J r , L | Lz ) , (K.5b)  \nwith F rot( J r , L | Lz ) = F rot( J r , L , Lz ) / F ( Lz ) standing for the DF in ( J r , L ) given Lz , and normalised to unity. Importantly, for Lz , 0, the coe ffi cients DL z and DL z Lz are well defined whatever the considered rotating DF. In addition, they both converge to some finite values as Lz → 0, as illustrated in Fig. K.3. The absence of any divergence here emphasizes that the divergence observed in Fig. K.2 stems from a coordinate singularity associated with cos I .', 'Appendix K.2: Impact on N -body measurements': 'To estimate the impact of the discontinuity on relaxation, we used N -body simulations in Fig. K.4 to compare the early relaxation of the discontinuous LBD distribution (equation 5) with its smooth approximation (equation K.2 with a = 10). Reassuringly, we observe a fast dilution of the discontinuity, with the DF quickly resembling its smooth approximation. Afterwards, no significant di ff erence is observed between the two systems.', 'K. Tep, J.-B. Fouvry & C. Pichon: Non-resonant relaxation': 'Fig. J.1. Same as in Fig. 9 but in ( L , cos I ). Di ff usion reshu ffl es orbital inclinations toward a more a ffi ne distribution in cos I .  \n<!-- image -->  \nFig. J.2. Same as Fig. 10 but in ( L , cos I ). The NR prediction fails to recover in detail the di ff usion structures observed numerically in Fig. J.1.  \n<!-- image -->  \nFig. K.4. Comparison between two rotating, isotropic Plummer clusters with parameters ( q , α ) = (0 , 0 . 25) and N = 10 5 stars. The discontinuous cluster follows the LBD parametrization (equation 6), while the smooth cluster follow an erf approximation, as defined in eq. (K.2). Both measurements have been ensemble-averaged over 10 realizations. Beyond the initial dilution of the discontinuity, we observe that the initially discontinuous DF joins the smooth DF, after which both systems appear to follow the same evolution.  \n<!-- image -->  \nFig. K.1. Family of functions ga (cos I ) (equation K.2) for various a . For a → 0, the function approaches identity, while for a →∞ , it approaches the sign function. Varying a allows us to investigate the impact of the LBD discontinuity in equation (5).  \n<!-- image -->  \nFig. K.2. Di ff usion flux along cos I of the 3D FP equation in ( J r , L , cos I ), for J r = 0 . 1 and cos I = 0 . 2, as a function of L and for various smooth functions ga (equation K.2). Here, we consider an isotropic cluster with rotation parameter α = 0 . 25. For any smooth ga , the flux diverges like 1 / L for L → 0. As a →∞ , that is, as ga tends to the sgn function, the flux converges to the LBD flux pointwise, which exhibits no divergence.  \n<!-- image -->  \nFig. K.3. One-dimensional di ff usion coe ffi cients DL z (top panel) and DL z Lz (bottom panel) of equations (K.5), computed for the rotating isotropic clusters with α = 0 . 25, and with the smoothing parameters a = 1 , 5 , 10 , 15 , ∞ , as defined in equation (K.2). These coe ffi cients are well-defined for all values of Lz . In addition, they only (very) weakly depend on a , with relative di ff erences of order 1-5%.  \n<!-- image -->'}