Nathan Magnan

Jean-Baptiste Fouvry

Christophe Pichon

Pierre-Henri Chavanis

Orbital alignment and mass segregation in galactic nuclei via vector resonant relaxation

Keywords: Diffusion -Gravitation -Galaxies, kinematics and dynamics -Galaxies, nuclei

Supermassive black holes dominate the gravitational potential in galactic nuclei. In these dense environments, stars follow nearly Keplerian orbits and see their orbital planes relax through the potential fluctuations generated by the stellar cluster itself. For typical astrophysical galactic nuclei, the most likely outcome of this vector resonant relaxation (VRR) is that the orbital planes of the most massive stars spontaneously self-align within a narrow disc. We present a maximum entropy method to systematically determine this long-term distribution of orientations and use it for a wide range of stellar orbital parameters and initial conditions. The heaviest stellar objects are found to live within a thin equatorial disk. The thickness of this disk depends on the stars' initial mass function, and on the geometry of the initial cluster. This work highlights a possible (indirect) novel method to constrain the distribution of intermediate mass black holes in galactic nuclei.

INTRODUCTION

Supermassive black holes (BHs) are ubiquitous in external galaxies (Magorrian et al. 1998;Genzel et al. 2010;Kormendy & Ho 2013) and their induced feedback plays a critical role in regulating galaxy formation through cosmic ages (Heckman & Best 2014). The unique proximity of our Galactic centre is an extraordinary opportunity to study and constrain the long-term evolution of galactic nuclei and the stellar clusters orbiting within. Recent developments in that realm include detailed census of stellar populations around SgrA* (Ghez et al. 2008;Gillessen et al. 2017), the observation of a cool accretion disc (Murchikova et al. 2019), as well as the observation of the relativistic precession of S2 (Gravity Collab. et al. 2020). Similarly, the origin of the observed clockwise stellar disc (Bartko et al. 2009;Yelda et al. 2014) has triggered a lot of interest, as its existence may impact the merger rate of the intermediate mass black holes (IMBHs) population [START_REF] Zwart | [END_REF].

As already pointed in [START_REF] Rauch | [END_REF], the steep potential well generated by the central BH allows for efficient secular orbit-averaged interactions between stars, driving their relaxation through an intricate hierarchy of dynamical processes (Alexander 2017). In this paper, we focus on the process of vector resonant relaxation (VRR) during which stars undergo a random reshuffling nathan.magnan@maths.cam.ac.uk of their orbital orientations through long-term coherent torques between the finite number of stellar orbital planes (Kocsis & Tremaine 2015). Given that VRR occurs on a timescale shorter than the stellar ages (see fig. 1 in [START_REF] Kocsis | [END_REF], one may expect that their observed distribution of orientations corresponds to some statistical equilibrium.

Determining the outcome of this long-term rearrangement has been the focus of recent efforts (Roupas et al. 2017;Takács & Kocsis 2018;Szölgyén & Kocsis 2018;Touma et al. 2019;Tremaine 2020a,b;Gruzinov et al. 2020) that jointly offer new clues on the fascinating properties of these long-range interacting systems, such as negative temperatures or phase transitions. In particular, Szölgyén & Kocsis (2018), using an ingenious Monte-Carlo approach, have shown that in systems with a wide range of stellar populations (i.e. various masses and semi-major axes), VRR can lead to the spontaneous formation of a disc through the angular segregation of the most massive stars and IMBHs. This is the issue that we further investigate in this work. We develop and implement an explicit and efficient optimisation procedure to find maximum entropy solutions compatible with an initial configuration, so as to infer the thermodynamical equilibria of a given galactic nucleus. In the particular context of VRR, this roadmap was already started out by Roupas et al. (2017) in the limit of a single-population system with a quadrupolar interaction, and later improved in Takács & Kocsis (2018) which, while still re-stricting themselves to a single-population system, considered harmonic expansion beyond the quadrupole. Here, we expand these works to multi-population clusters by emphasising the critical role played by the diversity of stellar orbits and masses to allow for nontrivial anisotropic statistical equilibria, as unveiled in Szölgyén & Kocsis (2018).

The paper is organised as follows. In §2, we briefly review the process of VRR in galactic nuclei. We also present the numerical optimisation method used to determine the thermodynamical equilibria. In §3, we use this approach to carefully investigate the details of equilibrium configurations that appear in systems composed of multiple stellar populations, as well as the impact of the initial mass function and the geometry of the formation scenario. Finally, we conclude in §4.

VRR MODEL

Let us first set up the framework in which to quantify the long-term effect of VRR on a population of orbital planes.

Interaction of Keplerian annuli via VRR

We consider an isolated cluster of N 1 stars orbiting a supermassive BH of mass M•. Following a double orbit-average over all the stars' Keplerian motion and in-plane precession, the VRR Hamiltonian takes the form1 (Kocsis & Tremaine 2015)

Htot = - N i<j Gmimj |ri(t) -rj(t )| t,t , (1)
where the sum over (i, j) runs over all the pairs of particles, and the double orbit-average • t,t operates over the fast Keplerian motions and in-plane precessions of both particles, with ri(t) and rj(t ) describing their trajectories.

Using the Legendre expansion of the Newtonian interaction, and the addition theorem for spherical harmonics, one can rewrite equation (1) as (see, e.g., Fouvry et al. 2019b)

Htot = - N i<j max =2 m=- H Ki, Kj Y m ( Li) Y m ( Lj), (2)
where Li stands for the unit vector aligned with the instantaneous orientation of the star's orbital plane and the real spherical harmonics, Y m ( L), are normalised so that d L Y m Y m = δ δ mm . In equation ( 2), a Keplerian annulus (see Fig. 1) is fully characterised by its conserved quantities

K = (m, a, e), (3)
with m the star's individual mass, a its semi-major axis, and e its eccentricity. We stress that, owing to orbit average, each Ki are taken to be fixed throughout the VRR dynamics. Finally, in equation (2), we introduced the (symmetric) coupling coefficient

H Ki, Kj = Gmimj 4π 2 +1 P (0) 2 × π 0 dMi π π 0 dMj π Min[ri, rj] Max[ri, rj] +1 , (4)
with P the Legendre polynomial of order and Mi the mean anomaly of the orbit Ki (see, e.g., [START_REF] Murray | Solar System Dynamics[END_REF]. These coupling coefficients fully encode the impact of the shape of the two orbits on their relative torque and show significant diversity (see, e.g., fig. 1

in Kocsis & Tremaine 2015) 2 .
Finally, we point out a few specificities of the multipole expansion from Eq. ( 2): (i) All odd harmonics have vanishing coupling coefficients, H = 0, so that they can be dropped; (ii) The Legendre expansion is truncated to the maximum order max;

(iii) We do not account for the harmonics = 0 as it does not drive any dynamics. With such a convention, spherically symmetric distributions have a vanishing mean total energy.

Describing the process of VRR amounts to describing the long-term evolution of Li, as driven by the Hamiltonian from equation (2) (see §2.2 in Kocsis & Tremaine 2015, for the associated equations of motion). We illustrate such a system in Fig. 1. One can equivalently describe the instantaneous state of the stellar cluster with its discrete distribution function (DF)

F d ( L, K, t) = N i=1 δD( L -Li(t)) δD(K -Ki), (5)
which follows the normalisation d LdKF d = N . Assuming that particles with the same orbital parameters are indistinguishable, this DF entirely describes the system's state. Therefore, characterising the VRR dynamics amounts to characterising ∂F d /∂t, using the Klimontovitch equation [START_REF] Klimontovich | The statistical theory of non-equilibrium processes in a plasma[END_REF]).

Equilibrium configurations

Fortunately, the VRR dynamics (∼ 1 Myr for S2 around SgrA*, see fig. 1 in [START_REF] Kocsis | [END_REF]) is rapid compared to the cluster's age (∼ 10 Myr for the S-cluster, see Habibi et al. (2017)).

As such, if one is interested in sufficiently long timescales, rather than describing the details of ∂F d /∂t, one may solely focus on characterising the expected equilibrium configurations reached at late times. This is given by

Feq( L, K) = lim t→+∞ F d ( L, K, t) , (6)
with • standing for an ensemble average over independent realisations of the system -see §2.3 for the appropriate handling of this ensemble average. Efficiently predicting these long-term equilibrium DFs is the focus of the present work. In such a late-time limit, the only information retained by the cluster are its invariants. These are:

• The number density of stars with orbital parameters K

N (K) = d L Feq( L, K). (7)
• The total energy

Etot = 1 2 d LdK Feq( L, K) ε( L, K), (8)
where ε( L, K) stands for the energy of a particle of orientation L and orbital parameter K as

ε( L, K) = - max =2 m=- dK H K, K M m (K )Y m ( L), (9)
with the magnetisations

M m (K) = d L Y m ( L) Feq( L, K) (10)
defined as spherical harmonic moments of the DF.

• The total angular momentum

Ltot = d L dK L(K) L Feq( L, K), (11)
with the norm of the angular momentum vector

L(K) = m GM• a (1 -e 2 ). (12)
As usually carried out in the microcanonical ensemble (i.e. for an isolated cluster), for a given set of invariants {N (K), Etot, Ltot}, the admissible equilibrium configurations are obtained by maximising the Boltzmann entropy

S = -kB d LdK Feq( L, K) ln Feq( L, K) , (13)
with kB the Boltzmann constant, under the previous conservation constraints. Such a maximisation generically yields (see §A)

Feq( L, K) = N (K) e -β ε( L,K)+L(K) γ• L d L e -β ε( L ,K)+L(K) γ• L , ( 14
)
where β and γ are the Lagrange multipliers respectively associated with the total energy and angular momentum conservation. We note that equation ( 14) is very similar to its co-planar counterpart found in the HMF model (Antoni & Ruffo 1995;Chavanis et al. 2005) to capture bar formation [START_REF] Pichon | Statistical Description of Transport in Plasma[END_REF] as a phase transition towards orbit alignment. Gruzinov et al. (2020) recently used a similar mean-field approximation and maximum entropy method to determine the thermodynamical equilibria of black hole star clusters. It is therefore no surprise that their equation ( 3) is so similar to the present equation (14). The main difference is that Gruzinov et al. (2020) considered the case of massive Keplerian elliptic wires, while we consider the case of massive Keplerian annuli. Phrased differently, Gruzinov et al. (2020) investigated the simultaneous equilibria of scalar resonant relaxation (SRR) and VRR -i.e. the joint relaxation of eccentricities and orientations -whereas we focus here on the equilibria of VRR -i.e. the sole relaxation of orientations.

Axisymmetric assumption

For a given nuclear cluster, the total angular momentum vector Ltot provides us with only one specific direction, taken to be the +z axis throughout the paper. Unfortunately, this does not imply that the thermodynamical equilibrium of VRR necessarily has an axial symmetry around Ltot, as a spontaneous symmetry breaking could occur (see, e.g., [START_REF] Kocsis | [END_REF]Gruzinov et al. 2020).

For simplicity however, we assume that VRR does not exhibit here any such symmetry breaking, and restrict ourselves to axisymmetric DFs, i.e.

∀ K, ∀ , ∀ m = 0, M m (K) = 0. ( 15
)
From there, we can make the simplifications Ltot → Ltot > 0 and γ → γ > 0. This assumption greatly reduces the total number of spherical harmonics to consider, therefore it significantly alleviates the numerical complexity. However, let us stress that this assumption is not always physically motivated. It is still legitimate in some cases, as for single-population and single-harmonic clusters (Roupas et al. 2017) or multi-population clusters with null inverse temperatures ( §C3), but it does not hold in some other regimes, as emphasised by the finding of a warped VRR disc in the simulations of [START_REF] Kocsis | [END_REF] (see fig. 6 therein).

Overall, the axisymmetric assumption is an important limitation of the present work, and needs to be challenged in future studies.

To comply with this approach, the ensemble average of equation (6) is carried out over realisations which are all rotated to have their Ltot aligned along +z.

Self-consistency

Of course, one needs to impose self-consistency on equations ( 8) and ( 11), as well as on equation ( 14). Indeed, Feq( L, K) involves the one-particle energy ε( L, K) which, via equations ( 9) and (10), involves Feq( L, K) itself. Within the present microcanonical ensemble, imposing the cluster's total energy, angular momentum and orbital distribution ultimately sets up its temperature (via β), rate of rotation (via γ), and shape (via M 0 (K))3 .

To effectively solve such a generic problem, we discretise the distribution as a finite set of stellar populations, indexed by k and described by the orbital parameters K k . We then write

Feq( L, K) = k F k ( L) δD(K-K k ) with F k ( L)
the distribution of orientation of the k th population and N k = d LF k ( L) its number of stars (see equation 7). After this discretisation, a cluster's configuration is fully characterised by its set of order parameters

θ = β, γ, {M ,k } , (16)
with M ,k = d LY 0 ( L)F k ( L) the axisymmetric -magnetisation of the k th population (see equation 10).

For a given initial condition, i.e. a given (Etot, Ltot, {N k }), the cluster's equilibrium configurations are obtained for the parameters θ that are joint roots of the consistency functions

CE = Etot -Etot Feq(θ) /Etot, CL = Ltot -Ltot Feq(θ) /Ltot, ∀ , ∀ k, CM ,k = M ,k -M ,k Feq(θ) /(N k y ), (17)
with y = (2 +1)/(4π) and Feq(θ) coming from equation ( 14). Importantly, in order to place all the constraints from equation ( 17) on equal footings, the consistency functions C = CE, CL, CM ,k are all dimensionless, and rescaled to be of order unity.

In order to find roots of the function θ → C[θ], we use Newton-Raphson's method (e.g., [START_REF] Press | Numerical Recipes 3rd Ed[END_REF]), following Takács & Kocsis (2018). More precisely, starting from a configuration θn, we compute

θn+1 = θn -J -1 [θn] C[θn], (18)
to obtain the next iteration4 , with J[θ]ij = ∂Ci/∂θj the Jacobian of the consistency function C at point θ. Given that C has a simple analytical form in equation ( 17), so does its Jacobian, as detailed in §B. Equation ( 18) is the second important difference between Gruzinov et al. (2020) and the present work. Indeed, Gruzinov et al. (2020) solves the self-consistency requirement on their equation ( 3) by iteratively computing a sequence ψ → F → ρ → ψ → • • •, of potentials, DFs and densities, up to convergence. Here, we directly ensure self-consistency in equation ( 14) by using Newton's method.

Optimisation strategy

To finalise our algorithm, it only remains to specify our choice for the starting point θ0 of the iteration process.

In §C1 we explore a first method where θ0 is initialised at random. Such an agnostic approach allows us to recover a cluster's both stable, metastable and unstable equilibria, and to extend the results of Takács & Kocsis (2018) to (axisymmetric) multipopulation clusters. Using this method, we recover all the qualitative behaviours of the single-population clusters reported in Roupas et al. (2017) and Takács & Kocsis (2018), and confirm that the branch that goes through β = 0 (see Fig. C1) always has the highest entropy. Unfortunately, as one increases the number of stellar populations, the efficiency of this protocol drastically drops.

Dealing with systems with many stellar populations therefore requires improvements to the initialisation process. To alleviate most of these difficulties, we restrict ourselves, and predict only the clusters' global thermodynamical equilibria, and none of the other possible equilibria, should they be unstable or metastable. Benefiting from the insight of Fig. C1, we obtain these global equilibria by iteratively moving along the series of equilibria associated with the branch that has a solution for β = 0 (see §C2). Of course, the main drawback of this approach is that it cannot be used to determine any of the unstable thermodynamical equilibria. and associated with distinct episodes of star formation or infall events. Each patch is made of Npart = 512 stars distributed uniformly within a small angular section around the patch's centre. Each star is coloured according to the norm of its angular momentum vector, with the smallest norm in blue and the largest one in red. Clearly, stars with small angular momenta are the most numerous.

THERMODYNAMICAL NUCLEAR EQUILIBRIA

We can now make use of our entropy optimisation algorithm to investigate the typical equilibrium distribution of orientations in galactic nuclei. Benefiting from the efficiency and versatility of this method, we will also use it to explore the (large) parameter space describing possible initial stellar clusters.

Parametrisation of the stellar population

Let us first parametrise the distribution of the orbital parameters K k , i.e., the distribution of masses, semi-major axes, and eccentricities. Our fiducial model is the same as in Szölgyén & Kocsis (2018). We assume that stars are formed through a series of distinct episodes of star formation or infall events. More precisely, for a given realisation, we generate N disc = 16 discs, each of them composed of Npart = 512 stars, so that the cluster's total number of particles is N = N disc ×Npart = 8 192. For each disc, its average orientation, L disc , is drawn uniformly on the the unit sphere, while the orientations of the stars within that disc are drawn uniformly within the small region L• L disc 0.994. Finally, for each star, the orbital parameters K = (m, a, e) are drawn independently from one another according to power law probability distribution functions (PDFs) proportional to (m -2 , a 0 , e) respectively within the ranges mmax/mmin = 100, amax/amin = 100, and (emin, emax) = (0, 0.3). In Fig. 2, we illustrate one typical realisation of such a protocol. Following §2.3, we recall that the realisations are always rotated to have their Ltot aligned along +z.

Once an initial distribution has been drawn, one may compute its two key invariants, Etot and |Ltot|. In practice, we keep track of these two invariants through the two dimensionless quantities

E = - Etot N 2 Gm 2 min /amin ; s = |Ltot| i L(Ki) , (19)
which we respectively call the binding energy and spin of the stellar cluster. In Fig. 3, we represent the typical distribution of (E, s) for a large number of stellar clusters drawn from our fiducial procedure. We find that clusters exhibit a wide range of binding energy and spin, and in turn this somewhat impacts the diversity of thermodynamical equilibria.

A typical equilibrium distribution

In Fig. 3, we determined the distribution of the clusters' invariants (E, s). In order to characterise the properties of the associated systems, we now investigate in detail the particular case (E, s) = (2×10 -3 , 0.2), i.e. a typical cluster realisation. We detail in §D1 our precise choices for the discretisation of the stellar populations. By maximising the entropy, we determine this cluster's equilibrium configuration, Feq( L, K). This is first illustrated in Fig. 4 in the (m, Lz)-plane. The same cluster is alternatively represented in Fig. 5, where we present the stellar density distribution on the sphere for two different mass bins.

As is clearly visible, heavy particles tend to have their unit angular momentum vector L oriented towards the North pole, meaning that the associated stars tend to orbit near the equatorial plane, a conclusion already reached by Szölgyén & Kocsis (2018) using an alternative Monte-Carlo approach5 , and by Gruzinov et al. (2020) in the more general case of Keplerian elliptic wires. Having the heavy particles, e.g., the IMBHs, orbit close to the same orbital plane should drastically impact the rate of their pairwise mergers in galactic nuclei. In practice, we also repeated the experiment from Fig. 4 a hundred times by letting (E, s) explore the distribution from Fig. 3. The associated ensemble-averaged DF was found to be similar to the one from Fig. 4.

For the heaviest stars, we also note the presence of an additional over-density near the South pole indicating the presence of a counter-rotating equatorial disc. However, owing to the conservation of angular momentum, this component is less populated than the main prograde disc. As one considers lighter particles, the anisotropy fades away and the lightest particles do not show any strong sign of spontaneous orientation alignment. This concurs with Gruzinov et al. (2020)'s results (see §6.3. therein), which found that light objects follow a spherically symmetric distribution.

Having the full equilibrium distribution Feq( L, K) at our disposal, one may also study its dependence w.r.t. the semi-major axis and eccentricity. Limiting ourselves to particles of intermediate mass, we find that the disc is slightly thinner for intermediate semi-major axes, as illustrated in Fig. 6. This is in agreement with

L  z F eq (a, L  z ) max L  z F eq (a, L  z )
a/a min = 1 a/a min = 10 a/a min = 100 As for eccentricities, we find that, in essence, Feq is independent of e. This was expected given that the coupling coefficients, H [K, K ] (equation 4) and the norm of the angular momentum L(K) (equation 12) only weakly depend on e, especially for the chosen quasi-circular orbits 0 e 0.3.

Impact of the cluster's binding energy and spin

As highlighted in Fig. 3, the clusters exhibit a significant diversity in their binding energies and spins, given our generation protocol. Let us therefore investigate the dependence of the equilibria w.r.t. these invariants. In Fig. 7, we present series of equilibria (caloric curves) giving the inverse temperature, β, as a function of the normalised total energy, -E, for various total angular momentum (i.e. various s). We observe that the temperature generically increases with energy, except in a small region of negative specific heat. Also, (m, a, e) drawn following the protocol of §3.1, and with binding energies distributed on a logarithmic grid ranging from 10 -6 to 10 -1 . Bottom panel: Illustration of the equilibrium angular frequency γ as a function of the cluster's spin for different binding energies. For each line, 100 systems were simulated, with spins distributed on a linear grid ranging from 0 to 1. See equation (C1) for the definitions of β 0 and γ 0 . As one lowers the total angular momentum, via s, thermodynamical equilibria can exhibit negative specific heats, i.e. C = ∂Etot/∂T < 0 (Roupas et al. 2017). For a given value of s, our approach to determine the series of equilibria is to start from the (much) simpler problem at β = 0 and move up to the target energy step by step iteratively (see §C2).

for any value of the spin there exists a region of negative temperature. Both effects had already been reported in Roupas et al. (2017) and Takács & Kocsis (2018), and are further discussed in §C1. Note however that these behaviours arise for E 10 -4 , i.e. well outside of the astrophysical regime, E 10 -3 , found in Fig. 3. In the same figure, we also present series of equilibria giving the angular frequency γ as a function of the total angular momentum s, for fixed values of the energy E. The angular frequency is found to always increase with the cluster's spin.

Within the domain of invariants spanned by Fig. 3, all equilibria are found to remain qualitatively similar but differ in the strength of their anisotropy. To characterise the level of anisotropy in the orbital distribution, we introduce a cluster's segregation rate, SR, as

SR = m mc dK Lz Lc d L Feq( L, K) m mc dK d L Feq( L, K) , (20)
This quantity describes which fraction of heavy particles (i.e. with m mc) have an angular momentum vector that is well-aligned with the cluster's total angular momentum (i.e. have Lz Lc). With such a definition, we note that anti-aligned particles are not accounted for. The larger SR, the stronger the anisotropic segregation of the heavy populations. In practice, guided by Fig. 4, we consider the values mc = 10 mmin and Lc = cos(20 • ). We expect that changing the values of mc does not qualitatively impact the present conclusions.

In Fig. 8, we illustrate the clusters' segregation rate as a function of the binding energy and spin (E, s). We observe the crucial role played by these two invariants in driving the spontaneous alignment of the orbital orientations of the heaviest particles. The larger E and the larger s, the larger the proportion of heavy particles that are aligned near the equatorial plane.

To strengthen this conclusion, let us finally determine the angular size, θc, of the northern polar cap that contains a fraction fc of the heavy particles' unit angular momentum vectors L. More precisely, we define θc through the implicit constraint

m mc dK Lz cos(θc) d L Feq( L, K) m mc dK d L Feq( L, K) = fc, (21)
noting once again that anti-aligned stars are not accounted for. Equivalently, θc is also the angular size of the equatorial disk that contains a fraction fc of the heavy particles. Using a fraction fc = 50%, the dependence of θc w.r.t. the clusters' invariants is illustrated in Fig. 8. Once again, the more bound the cluster, and the larger its total angular momentum, the thinner the disc of heavy particles in the relaxed equilibrium.

Interestingly, if equations ( 20) and ( 21) were to consider both aligned and anti-aligned stars, through the replacement Lz cos(θc) → | Lz| cos(θc), the dependence w.r.t. s in Fig. 8 would be significantly reduced. This result is expected since in equation ( 8), the total energy is left invariant by the changes Li → -Li, while in equation ( 11) the total angular momentum changes sign.

Impact of the mass and semi-major axes distributions

Let us now step out and explore the impact of the distribution of the orbital stellar parameters themselves on the clusters' thermal equilibria. Following §3.1, we now vary the PDFs w.r.t. which the stars' masses and semi-major axes are drawn, keeping all other parameters the same. More precisely, we still draw (m, a) pairs according to power law distributions, (m γm , a γa ), but this time vary the power law indices (γm, γa) between different realisations.

The impact of changing these two parameters on the average segregation rate of the heavy particles, θc , is illustrated in Fig. 9. We find that the slope γm has the strongest effect with clusters containing fewer heavy stars segregating in a disk twice thinner than the fiducial cluster. The semi-major axis power law index also has an impact on the segregation angle: clusters with larger semimajor axes tend to display a stronger mass segregation. These dependencies may possibly be used in the future to estimate indirectly the stellar and compact objects' initial mass function in galactic centres. Indeed, the rate of pairwise IMBH mergers should be measurable using gravitational waves, and could serve as an indirect probe for the width of the segregated massive disc, which we just have shown is linked to the IMBHs' initial mass function.

Impact of the dispersion of disc orientations

A cluster's long-term equilibrium distribution is fully characterised by its three invariants: the binding energy, E, the spin, s, and the distribution of orbital parameters, N (K). Unfortunately, while dynamically relevant, such parameters do not translate easily in terms of astrophysical observables. As such, let us finally switch to Here, for each value of (γa, γm), the average has been computed over 128 realisations of a cluster. The fiducial case from Fig. 4 corresponds to (γa, γm) = (0, -2). Interestingly, both power law indices are found to impact the relaxed segregation angle.

parametrisations more closely related to the underlying formation process of the galactic centre's stellar cluster, such as the geometry of the initial distribution. Building upon Fig. 2, rather than drawing the initial orientations of the discs uniformly on the unit sphere, we may draw them along some biased direction, for example to reflect the preferential infall of new star-forming gas along specific directions, e.g., imposed by the geometry of past gaseous accretion events. More precisely, we assume that the average orientation of each disc is drawn according to a Von Mises-Fisher PDF (Wood 1994) of the form

P ( L) = κ 4π sinh(κ) e κ Lz , (22)
with κ the PDF's concentration. Here, κ = 0 corresponds to the isotropic case considered in Fig. 2. The larger κ, the smaller the spread of the PDF on the unit sphere, as in Fig. 10 for κ = 5. Keeping all other stellar parameters as in §3.1, we may now study the impact of the concentration parameter κ on the shape of the relaxed distributions. This is first illustrated in Fig. 11, where we show the dependence of the invariants (E, s) with κ. As expected, we recover that the more concentrated the distribution of the discs' initial orientations, the larger the value of s and E. In Fig. 12, we subsequently illustrate the dependence of the average segregation angle, θc (equation 21), as a function of κ. Narrower initial distributions of orientations lead to stronger segregation of the heavy stars at equilibrium. This is expected, since Fig. 8 showed that the segregation strength correlates positively with both binding energy and spin. To emphasise this conclusion, Fig. 12 also illustrates the thickness of the asymptotic massive disc as a function of the initial mass function slope. Similarly to Fig. 9, at fixed initial anisotropy, the steeper the slope, the stronger the mass segregation. Overall, investigations as in Fig. 12 should ultimately prove useful to place some constraints on the origin of SgrA*'s surrounding stel- . Typical distribution of the spin and binding energies as one varies the concentration κ of the initial discs. Coloured points correspond to the mean values of (s, E), while black error bars are standard deviations in both parameters. Statistics are computed from 10 4 realisations. As expected, the more concentrated the discs' orientations, the more packed the initial distributions, and therefore the larger the spins and binding energies. lar cluster, in particular in the light of the "paradox of youth" (Ghez et al. 2003;Genzel et al. 2010).

CONCLUSION

In the spirit of Roupas et al. (2017); Takács & Kocsis (2018), we relied on maximum entropy methods to assess the internal stellar structure of galactic nuclei and their underlying distribution of orbital orientations. We expanded their approach to multi-populations clusters and showed how to jointly account for the constraints of energy, angular momentum, and orbital parameters conservations to efficiently characterise a cluster's expected thermodynamical equi-

0 1 2 3 4 5 0 30 60 90 κ 〈θ c 〉 [°] γ m = -2.5 γ m = -2.0 γ m = -1.5
Figure 12. Average segregation angle, θc , as a function of the concentration parameter κ, for three different mass distributions. The solid lines represent the average segregation angle, and the dashed lines represent the standard deviation. The statistics are estimated from 2 10 realisations -though ∼ 10% of the realisations for (κ, γm) = (5, -2.5) failed because of numerical overflows encountered with very strongly anisotropic DFs. The more anisotropic a cluster initially, the higher its binding energy and spin, and therefore the more segregated at equilibrium. libria. In particular, as already pointed out in Szölgyén & Kocsis (2018), we recovered the spontaneous alignment of the heavier stellar components, such as IMBHs (Fig. 4). The more bound the cluster and the larger its spin, the thinner the asymptotic disc (Fig. 8).

Benefiting from the versatility of the present method, we used it to explore the dependence of these structures w.r.t. the cluster's initial conditions of formation, such as the intrinsic power law spectra in mass and semi-major axis (Fig. 9), or the initial spread in orientations (Fig. 12). Small modifications to the present algorithm would also make it possible to confirm that the final equilibrium does not depend on the specific sequence of star formation events (Szölgyén & Kocsis 2018).

Of course, this work is but one step towards the systematic exploration of the remaining imprints of a galactic nucleus' initial formation on the cluster's long-term distribution of orientations. Let us conclude by listing a few possible venues for future works.

The present approach could first benefit from various extensions on the theoretical and analytical front. (i) We have assumed from the start that stellar eccentricities were conserved throughout the VRR evolution. This amounts to neglecting the process of SRR [START_REF] Rauch | [END_REF], whose signature is clearly visible in the S-cluster distribution (see, e.g., Tep et al. 2021). As already recently hinted in Szölgyén et al. (2021), allowing for the eccentricities to also vary during the orientations' relaxation might impact the system's VRR equilibria. This deserves careful study. Similarly, stars can also change their semi-major axis through non-resonant relaxation (NR), whose impact on the outcome of VRR should also be investigated. (ii) We recovered here that VRR may lead to highly non-spherical distribution of orientations. Such an end-state distribution is expected to impact in turn the rate of orientation diffusion of a given test star and the efficiency with which newly formed stellar discs can dilute (see, e.g., Giral Martínez et al. 2020).

The stochastic process of spontaneous alignment is not limited to the axisymmetric relaxation of Keplerian annuli around supermassive BHs. As such, at a significantly larger numerical cost, one could extend the present approach to non-axisymmetric distributions, as was already done in Roupas et al. (2017) in the limit of a single-population system with a quadrupolar interaction. Similarly, it would be of interest to investigate the distinction between the present global equilibria and the other possible metastable ones6 : do they also exhibit a disc-like structure? Finally, one should investigate whether there exist astrophysically relevant regimes with negative temperatures or heat capacities (Figs. 7 andC1).

In equation ( 1), we assumed that we could average the cluster's Hamiltonian over the in-plane precessions. We could lift this assumption and determine the equilibrium distributions of Keplerian elliptic wires (see, e.g., Gruzinov et al. 2020;Tremaine 2020b).

Finally, since the mean potential remains on average spherically symmetric, all these investigations could also be performed in globular clusters, where stellar metallicities or ages could be used as additional tracers of the mixing of stellar orientations.

In the present work, we restricted ourselves to predicting the end-state of VRR, i.e. its thermodynamical equilibrium. By design such an approach cannot provide any estimate of the expected relaxation time required for this asymptotic distribution to be reached. There are, at least, two possible venues to quantify such relaxation time: (i) One could use direct time integrations of the equations of motion to get the equilibrium distributions, a goal already pursued in Kocsis & Tremaine (2015). Building upon Fouvry et al. (2020), one can expect that efficient multipole methods may be designed to perform such direct numerical simulations with a computational complexity scaling linearly with the total number of particles; (ii) In the limit of sufficiently symmetric orbital distributions, e.g., axisymmetric (Fouvry et al. 2019a), one could alternatively derive an explicit kinetic theory for the cluster, from which the relaxation time would naturally follow. Interestingly, we report that the axisymmetric equilibrium recovered in Fig. 4 generically exhibits a monotonic profile of latitudinal precession frequency for all stellar populations. This may induce a situation of "kinetic blocking" and play a critical role in defining the efficiency with which these systems may relax (Fouvry et al. 2019a).

Finally, benefiting from the planned upgrade on VLTI [START_REF] Eisenhauer | The Very Large Telescope[END_REF]Gravity Collab. et al. 2021), as well as the future thirtymeter class telescopes such as ELT [START_REF] Pott | Ground-based and Airborne Instrumentation for Astronomy VII[END_REF][START_REF] Davies | Ground-based and Airborne Instrumentation for Astronomy VII[END_REF] and TMT [START_REF] Do | [END_REF]), we will soon have a wealth of orbitally-resolved stars around SgrA* along with their stellar ages. Building upon the present work, a detailed characterisation of their orbital distribution should prove paramount to place constraints on the properties of the (likely present) un-observed IMBHs, and the impact of their effective distribution of orientations on their overall in-situ merger rates.

APPENDIX A: THERMODYNAMICAL EQUILIBRIUM Following Roupas et al. (2017), we derive the generic expression of the DF of thermodynamical equilibria, Feq, from equation ( 14). This DF maximises the Boltzmann entropy (equation 13), under three joint constraints: the conservation of orbital parameters (equation 7), the conservation of energy (equation 8) and the conservation of angular momentum (equation 11).

Following the double orbit-average from equation ( 1), we stress that the Hamiltonian of VRR (equation 2) is such that (i) the phase space domain for L is of finite volume and (ii) for a finite max, there is no divergence of the pairwise VRR interaction for Li = Lj. As such, the VRR dynamics is simplified and does not present the typical peculiarities of the statistical mechanics of self-gravitating systems such as the evaporation of stars, the formation of binaries, the gravothermal instability, and the absence of a true statistical equilibrium (see, e.g., Chavanis 2006). In the present case, we essentially deal with the dynamics of self-interacting spins on a sphere, a system which possesses well-defined statistical equilibria in all cases.

Let us consider Feq a local extremum of entropy. Lagrange multipliers then guarantee the existence of a function α(K), a scalar β and a vector γ that ensure the differential equality DF eq S + dK α(K)DF eq N (K)-βDF eq Etot +γ •DF eq Ltot = 0.

(A1)

In that expression, the differentials are given by the linear forms

DF eq S : δF → -kB d LdK 1+ln Feq( L, K) δF ( L, K), DF eq N (K) : δF → d L δF ( L, K), DF eq Etot : δF → d LdK ε( L, K) δF ( L, K), DF eq Ltot : δF → d LdK L(K) L δF ( L, K), (A2)
where we used the symmetry

H [K, K ] = H [K , K ] (see equa- tion 4
) to compute DF eq Etot. Injecting these relations into equation (A1), we find

∀ δF, d LdK -kB 1 + ln Feq( L, K) + α(K) (A3) -β ε( L, K) + L(K) γ • L δF ( L, K) = 0,
Since this integral must vanish whatever the small displacement δF considered, we have

∀ L, ∀ K, -kB 1 + ln Feq( L, K) + α(K) + β ε( L, K) + L(K) γ • L = 0. (A4)
Inverting this relation, we find that the equilibrium DF is necessarily of the form

Feq( L, K) = exp α(K)-1 -β ε( L, K)+L(K) γ • L , (A5)
where, for convenience, the factor kB has been absorbed in the definition of α, β and γ. Following equation ( 7), we impose the number density of stars with orbital parameters K to be N (K), and we finally recover equation ( 14). In Fig. C1, we use this protocol for a multi-population and a multi-harmonic cluster with (Nm, Na, Ne) = (5, 1, 1), and max = 10 -see §D for details on the cluster's properties. In a sense, this figure is similar to the top panel of Fig. 7, except that in Fig. C1, we do not restrict ourselves solely to the global thermodynamical equilibria. Interestingly, it also appears that all the qualitative behaviours reported in Roupas et al. (2017) and Takács & Kocsis (2018) are still present in the multi-population case, as can be seen by comparing with Figs. D3 andD4.

We first note that for any spin s, there exists a region of the caloric curve where equilibria have negative temperatures. Negative temperatures occur when the total volume of phase space is finite (which is the case for VRR), as first shown by [START_REF] Onsager | [END_REF] in the context of 2D vortex dynamics. However, in the present context, the negative temperature states do not seem to have any particular impact and, furthermore, the region β < 0 is found to be outside of the astrophysically relevant regime (see Fig. 3).

For some binding energies and spins, there can exist several solutions to equation ( 14). Indeed, the Lagrange multipliers method finds all local extrema of the entropy. These could be local minima, saddle points, metastable local maxima, or the one global maximum. In the present context, the turning point method of Poincaré (see, e.g., Lynden-Bell & Wood 1968;Katz 1978) states that along a continuous branch in the (Etot, β)-plane, there can be a loss of stability around a given solution in the microcanonical ensemble (i.e. fixed Etot and Ltot) only if Etot → β(Etot) exhibits an infinite derivative there. In Fig. C1, no such turning points are observed. Given that the solution for β = 0 is always found to be stable (see §C3), we can conclude that the branch of solutions that goes through β = 0 are always, at least, metastable. In practice, we systematically computed the entropy of the various branches, and the one going through β = 0 was always found to correspond to the global (axisymmetric) maxima. This justifies our choice in §C2 to limit ourselves to only determining the equilibria along β = 0, therefore significantly alleviating the numerical burden.

In Fig. C1, we also recover that low-spin multi-population clusters can exhibit negative specific heats at small energies. However, when considering invariants (E, s) provided by the protocol from §3.1, none of the realisations considered were found to fall within that domain. Finally, we note that the region of negative specific heat is delimited by a pair of equilibrium points where the function Etot → β(Etot) has a vanishing derivative. The Poincaré turning point method can be used to show that the series of equilibria in between these two points of vanishing derivative is unstable in the canonical ensemble (fixed β and γ). This leads to a situation of ensemble inequivalence -since equilibria with negative specific heats are allowed in the microcanonical ensemble but forbidden in the canonical ensemble -characteristic of systems driven by longrange interactions (see, e.g., Chavanis 2006). Given the similarities of Fig. C1 with the results already presented in Roupas et al. (2017); Takács & Kocsis (2018), we refer to these previous works for detailed discussions of these canonical phase transitions. We stress, however, that the relevant ensemble considered in the main text is the microcanonical one since our system is assumed to be isolated.

C2 Iterative resolutions

When dealing with systems with a much larger number of stellar populations, the random initialisation from equation ( C1) is not efficient anymore. Indeed, because the total number of order parameters scale like O(Npop max) (see equation 16), the likelihood of starting sufficiently close to one equilibrium for equation ( 18) to converge drastically drops. As such, benefiting from the insight of Fig. C1, we limit ourselves to predicting the cluster's global thermodynamical equilibrium, i.e. predicting the equilibrium that lies on the branch that goes through β = 0.

More precisely, assuming that the cluster at hand is characterised by the two invariants (Etot, Ltot), we first determine the energy E (0) tot and order parameters θ (0) associated with the thermodynamical equilibrium at temperature β = 0 and total angular momentum Ltot. We highlight in §C3 how such a problem is significantly easier to solve, as it does not involve any self-gravitating contribution, therefore making equation ( 14) explicit.

From this initial configuration, we use θ (0) as an initial condition to solve a new self-consistency problem for the invariants (E tot the Newton method from §2.4 is expected to converge rapidly and to provide us with a new configuration θ (1) , solution of the problem with invariants (E (1) tot , Ltot). This new solution may then be used as an appropriate initial condition to solve the problem (E (2) tot , Ltot), getting us closer to our target energy Etot. Repeating this iterative process, we can ultimately solve the problem truly at hand, i.e. the one associated with the invariants (Etot, Ltot). In practice, we generically used a total of 50 different energies E (i) tot spread logarithmically between E (0) tot and the target total energy Etot.

C3 Equilibria for β = 0

We derive the DF that solves the entropy maximisation problem, without any constraint on energy. This DF is the starting point of our optimisation strategy with small steps in total energy (see §C2).

First, for a given non-zero value of Ltot (taken to be along +z), the solution of the optimisation problem without energy constraint must be axisymmetric. Indeed, if F (φ, Lz, K) satisfies the constraints from equations ( 7) and ( 11), then so does F∆ : (φ, Lz, K) → F (φ+∆, Lz, K). Given that equation ( 11) is linear w.r.t. F ,

F : (φ, Lz, K) → 2π 0 d∆ 2π F∆(φ, Lz, K) = F ( Lz, K) (C2)
also meets the constraints on Ltot. We note that F is a barycentre of the phase-shifted DFs, which all have the same entropy (equation 13). Given that s : x → x ln(x) is a convex function, Jensen's inequality gives S[F ] S(F ). As such, in the absence of any constraint on energy, the entropy maximum must be axisymmetric. Given this axisymmetry and remembering that β = 0 since there is no constraint on the energy, equation ( 14) becomes

Feq( L, K) = N (K) exp L(K) γ Lz d L exp L(K) γ L z . (C3)
Equation ( C3) is much simpler than equation ( 14) because its r.h.s. does not contain Feq anymore, i.e. the equation becomes explicit.

To find a value of γ for which Feq meets the constraint on Ltot, we compute the integral on the r.h.s. of equation ( 11). We get

L(K) d L Lz exp L(K) γ Lz = 2π L(K) +1 -1 dz z e L(K)γ z = 4π γ cosh[L(K) γ] -sinch[L(K) γ] , (C4)
and

d L exp L(K) γ Lz = 2π +1 -1 dz e L(K)γ z = 4π sinch[L(K) γ], (C5)
with sinch(x) = sinh(x)/x. Therefore, equation (11) finally gives

Ltot(γ) = dK N (K) L(K) cotanh[L(K) γ] -(1/γ) . (C6)
This is a strictly increasing function of γ, so that there exists at most one solution γ for a given value of Ltot. Although not analytical, such a solution is straightforward to obtain by dichotomy.

APPENDIX D: NUMERICAL APPLICATIONS

We briefly detail some of our choices in the effective numerical implementation of the method of entropy optimisation from §2.4 We also validate it in single-population systems by reproducing previous published results.

D1 Discretisation of the stellar populations

As explained in §2.4, in order to be effectively implemented the entropy optimisation requires a discretisation of the stellar populations. Each population is associated with an index k and is characterised by some orbital parameters K k and a particle number N k . In §3.1, we assumed that (m, a, e) are drawn independently from one another, so that they may be discretised independently as well. For the stellar mass, the considered range mmin m mmax is discretised in Nm logarithmic bins. Similarly, for a we use Na logarithmic bins, and Ne linear bins for e. As a consequence, the effective number of populations is set by Npop = NmNaNe. Then, for a given population, the number of particles N k is simply set by the value of the sampling power-law PDFs in the centre of the bins, multiplied by the volumes of the bins. Similarly, the value of the . Illustration of the dependence of (E, s) as a function of max for the protocol from §3.1. For a given value of max, the two plotted contours correspond to 33% (dashed) and 66% (solid) of the PDF's maximum. As already pointed out in Takács & Kocsis (2018), for max 10, the values of the two invariants (E, s) can be considered as converged.

population's orbital parameter, K k , is set by the value of K in the centre of the bin. As such, all our integrals over dK are formally approximated by Riemann sums using the midpoint rule.

To obtain Figs. 4 and 5, we used (Nm, Na, Ne) = (20, 10, 5). For the subsequent sections, we considered the simpler regime (Nm, Na, Ne) = (10, 10, 3). We checked that such choices of population numbers did not affect the results. Given that the numerical complexity scales like O(N 3 pop ), this greatly eases the overall parameter exploration.

D2 Convergence w.r.t. max

In practice, the Hamiltonian from equation (2) has to be truncated to some finite harmonic order max. While Roupas et al. (2017) limited themselves to the quadrupolar case max = 2, Takács & Kocsis (2018) lifted this restriction in single-population systems and showed that an effective truncation at max = 10 is sufficient to get converged results. This is what we briefly explore for multipopulation clusters in this section.

First, in Fig. D1, considering the same initialisation protocol as in §3.1, we illustrate the dependence of the invariants (E, s) as a function of max. In that figure, we recover that the typical value of E increases with max, while, of course, the average value of s is independent of it. In addition, we note that restricting oneself to max = 2 seems insufficient, while max = 10 offers reasonably well converged values of the binding energy.

In addition to affecting the values of the cluster's invariants, increasing the value of max might also impact the overall shape of the reconstructed anisotropic equilibria DF. This is what we explore in Fig. D2, where we compute the relative error in the reconstructed DF between max = 10 and max = 50. Given that the maximum relative error is ∼ 7% for the present clusters, we systematically truncated the pairwise interaction at max = 10 in all the figures presented in the main text. As the complexity of the en- tropy optimisation scales like O( 3 max ), this significantly alleviates the numerical difficulty of the computations.

D3 Validation in single-population systems

We validate our implementation of entropy maximisation by recovering previous results from the literature.

First, in Fig. D3, we recover the caloric curve presented in fig. 11 of Roupas et al. (2017) for a single-population system interacting only through = 2. In order to exactly match our present Takács & Kocsis (2018). In addition to the presence of negative specific heats for low s, we recover for weakly bound clusters the presence of more than one equilibrium solutions. Configurations with the smallest β were always found to have the highest entropy, hence correspond to the cluster's (axisymmetric) global statistical equilibria.

normalisation convention with the ones from Roupas et al. (2017), one has to consider βRoupas = β 3 8 N Gm 2 /a; ERoupas = Etot/ 3 8 N 2 Gm 2 /a . (D1) In this figure, it is interesting to recall that clusters with sufficiently small total angular momentum can present a negative specific heat.

In Fig. D4, we recover fig. 4 from Takács & Kocsis (2018), for a single-population system coupled beyond the quadrupolar interaction. To exactly match the normalisation conventions from Takács & Kocsis (2018), one has to consider β Takács = β N Gm 2 /a; E Takács = Etot/ N 2 Gm 2 /a . (D2)

In addition, following equation (2) of Takács & Kocsis (2018), we also had to replace the coupling coefficients from equation (4) with the simpler asymptotic scaling

H Takács = Gm 2 a 4π 2 (2 + 1) . (D3)
For sufficiently small values of the total spin s, one recovers more than one equilibrium solutions. The solutions with the smallest β, i.e. the branch that goes through β = 0, were always found to have the largest entropy, hence corresponding to the (axisymmetric) thermodynamical equilibria.

Figure 1 .

 1 Figure 1. Illustration of the process of VRR in a galactic nucleus. Following an orbit-average over the fast Keplerian motion induced by the central BH and the in-plane precessions (due to both the stellar mean potential and the relativistic corrections), stars are replaced by massive annuli (here quasicircular) which can torque one another (see also fig. 1 in Giral Martínez et al. 2020). This leads ultimately to the relaxation of the stellar orbital orientations with possibly the spontaneous formation of aligned discs.

Figure 2 .

 2 Figure2. Illustration of a cluster's typical initial condition. It consists of N disc = 16 dense patches of stars scattered uniformly on the unit sphere, and associated with distinct episodes of star formation or infall events. Each patch is made of Npart = 512 stars distributed uniformly within a small angular section around the patch's centre. Each star is coloured according to the norm of its angular momentum vector, with the smallest norm in blue and the largest one in red. Clearly, stars with small angular momenta are the most numerous.

Figure 3 .

 3 Figure 3. Distribution of the spin and binding energy of clusters drawn following the fiducial protocol from §3.1, for 2 18 realisations of initial conditions and max = 10.

FFigure 4 .

 4 Figure4. Top panel: Illustration of the relaxed stellar density distribution, Feq(m, Lz) (integrated over all semi-major axes and eccentricities), as a function of mass and orientation for a relaxed cluster with (E, s) = (2×10 -3 , 0.2). See §D for the numerical details. Our normalisation allows for each mass bin to have its maximum equal to 1. Heavier populations are much more segregated towards the poles. Bottom panel: Same as the top panel, where each line corresponds to a different mass bin.

Figure 5 .

 5 Figure 5. Stellar density distribution on the unit sphere, after relaxation, for the lightest particles (left, 1 m/m min 1.25) and the heaviest ones (right, 79 m/m min 100).We use the same normalisation and colouring as in Fig.4. Importantly, we note that the heavy particles are significantly more segregated towards the poles, i.e. the corresponding heavy stars live within a thin equatorial disc.

Figure 6 .

 6 Figure6. Illustration of the relaxed stellar density distribution, Feq(a, Lz) (integrated over all eccentricities, but only intermediate masses in 2 m/m min 16), as a function of semi-major axis and orientation, for the same cluster as in Fig.4. The normalisation allows for each mass bin to have its maximum equal to 1. Populations with an intermediate semi-major axis tend to segregate within a thinner disk.

Figure 7 .

 7 Figure 7. Top panel: Illustration of the equilibrium inverse temperature β as a function of the cluster's normalised total energy, -E (equation19) for different values of spin. Each line is composed of 500 systems, all with (m, a, e) drawn following the protocol of §3.1, and with binding energies distributed on a logarithmic grid ranging from 10 -6 to 10 -1 . Bottom panel: Illustration of the equilibrium angular frequency γ as a function of the cluster's spin for different binding energies. For each line, 100 systems were simulated, with spins distributed on a linear grid ranging from 0 to 1. See equation (C1) for the definitions of β 0 and γ 0 . As one lowers the total angular momentum, via s, thermodynamical equilibria can exhibit negative specific heats, i.e. C = ∂Etot/∂T < 0(Roupas et al. 2017). For a given value of s, our approach to determine the series of equilibria is to start from the (much) simpler problem at β = 0 and move up to the target energy step by step iteratively (see §C2).

Figure 8 .

 8 Figure8. Top panel: Segregation rate of the relaxed clusters (equation 20), within the range of spins and binding energies expected for typical astrophysical clusters (see Fig.3). Bottom panel: Same domain as the top panel but for the segregation angle θc (equation 21). Clusters with large binding energies and spins exhibit the strongest segregation of their heavy particles within the equatorial disc.

Figure 9 .

 9 Figure 9. Average segregation angle, θc , as a function of the power law indices of the initial distributions in semi-major axes and masses.Here, for each value of (γa, γm), the average has been computed over 128 realisations of a cluster. The fiducial case from Fig.4corresponds to (γa, γm) = (0, -2). Interestingly, both power law indices are found to impact the relaxed segregation angle.

Figure 10 .

 10 Figure 10. Same as Fig. 2, but assuming that the orientations of the discs are drawn from a von Mises-Fisher PDF with concentration κ = 5 (equation 22).

Figure C1 .

 C1 Figure C1. Inverse temperature β as a function of the normalised total energy, -E (see equation 19), for various spins s, when considering a cluster composed of 5 populations of different individual masses interacting with max = 10. This figure was obtained without restricting ourselves to the global thermodynamical equilibria (see C1), as emphasised by the second set of (red and blue) solutions at high β.

  Figure D1. Illustration of the dependence of (E, s) as a function of max for the protocol from §3.1. For a given value of max, the two plotted contours correspond to 33% (dashed) and 66% (solid) of the PDF's maximum. As already pointed out inTakács & Kocsis (2018), for max 10, the values of the two invariants (E, s) can be considered as converged.

Figure D2 .

 D2 Figure D2. Relative error in Fig. 4 as one increases the harmonic truncation from max = 10 to max = 50. The maximum relative error is ∼ 7%.

Figure D3 .

 D3 Figure D3. Inverse temperature β as a function of the energy E Roupas , for various fixed total angular momentum, s, for a single-population cluster with max = 2. This figure reproduces fig. 11 of Roupas et al. (2017). For low s, the cluster can exhibit a negative specific heat, i.e. ∂β Roupas /∂E Roupas > 0.

Figure D4 .

 D4 Figure D4. Same as Fig. D3 but considering max = 58. This figure reproduces fig. 4 ofTakács & Kocsis (2018). In addition to the presence of negative specific heats for low s, we recover for weakly bound clusters the presence of more than one equilibrium solutions. Configurations with the smallest β were always found to have the highest entropy, hence correspond to the cluster's (axisymmetric) global statistical equilibria.

We also assume that N i=1 m i M• (Roupas

2020).

Following §A ofFouvry et al. (2019b), the present coupling coefficients simply read H [K, K ] = J [K, K ] L(K), with L(K) introduced in equation (12). In practice, following the notations from §A ofFouvry et al. (2019b), the coupling coefficients were pre-computed once, up to max = 50 and on a 200×100×100 linear grid in (ln(a in /aout), e in , eout) with 10 -2 a in /aout 1 and 0 e in , eout 0.99, to be subsequently interpolated.MNRAS 000, 000-000 (0000)

In principle, there are infinitely many harmonics , but in practice we found it acceptable to stop at max = 10, see §D2.MNRAS 000, 000-000 (0000)

Note that in practice, it is faster and numerically more stable to solve the linear equation J x = -C for the unknown x = θ n+1 -θn, rather than to explicitly compute the inverse matrix, J -1 .

MNRAS 000, 000-000 (0000)

Note that even when accounting for the different normalisation, the details of the DF in (m, Lz)-space obtained here in Fig.4somewhat differ from the equivalent fig.2ofSzölgyén & Kocsis (2018).

For systems with long-range interactions, metastable states have a very long lifetime scaling as e N (see, e.g.,Chavanis 2006), because a system trapped in a metastable state (local entropy maximum) has to cross a huge barrier of potential to reach the fully stable equilibrium (global entropy maximum).MNRAS 000, 000-000 (0000)

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é. NM acknowledges support from SU-PAERO, the Université Paul Sabatier and the École polytechnique. We thank Stéphane Rouberol for running smoothly the Infinity cluster, where the simulations were performed.

DATA AVAILABILITY

The data and numerical codes underlying this article were produced by the authors. They will be shared on reasonable request to the corresponding author. The code is distributed on Github at the following URL: https://github.com/NathanMagnan/VRROOMe.

APPENDIX B: CONSISTENCY FUNCTIONS AND JACOBIAN

We detail here the expressions of the consistency functions from equation (17). In order to shorten the notations, let us first write the exponential factor from equation ( 14) as

We may then introduce the four angular integrals

The consistency functions from equation ( 17) then simply read

As required by equation ( 18), it is straightforward to compute the gradients of (CE, CL, {CM ,k }) w.r.t. the order parameters (β, γ, {M ,k }). All these expressions are analytical and involve angular integrals similar to the ones present in equation (B2). In practice, having discretised the stellar populations (see §D1), the integrals w.r.t. dK simply become k . In addition, benefiting from the assumption of axisymmetry (see §2.3), the integrals w.r.t. d L become integrals w.r.t. d Lz, which are computed using a Gauss-Legendre quadrature (see, e.g., [START_REF] Press | Numerical Recipes 3rd Ed[END_REF]) with 100 nodes.

APPENDIX C: OPTIMISATION STRATEGIES

In this Appendix, we present our two main approaches to solving the iteration problem from equation ( 18).

C1 Random initialisations

A first approach to initialise equation ( 18) is to consider starting points θ0 taken at random. More precisely, we draw the initial order parameters independently from one another and uniformly within the domains

with y = (2 +1)/(4π). The main advantage of such an agnostic approach is that it allows us to recover a cluster's both stable and unstable equilibria.


N.-A Masarra

J.-C Quantin

M Batistella

R El Hage

M F Pucci

R Ravel

J.-M Lopez-Cuesta

evolution of the volume electrical conductivity as a function of the PLA content A B

Introduction

• Percolation theory: Explains the transition from an insulating material to a conductive one in conductive filled polymer composites.

• Electrical percolation threshold: Minimal filler concentration required to acquire a conductive network within the polymer matrix • How to reduce the electrical percolation threshold? By creating a double electrical percolation threshold resulting from adding the conductive filler to an immiscible polymer blend matrix Required conditions for the double electrical percolation threshold

Compression molded samples

Selective localization of graphene in the PCL phase Co-continuous microstructure

Processing steps:

• Twin-screw extrusion to blend PLA, PCL, and GNP together

• Compression molding to prepare samples for analysis

• PLAx/PCLy/10 wt.% GNP were produced (x varies between 0 wt.% and 100 wt.% in the polymer total weight percentage). • PLA/PCL/GNP compression molded composites: Permanent selective localization of GNP in the PCL phase (The graphene + PCL zones are indicated by yellow arrows.)

• Good dispersion of graphene particles in the PCL phase

• Co-continuous microstructure is present in PLA65/PCL35/10 wt.% GNP and PLA60/PCL40/10 wt.% GNP formulations (The brighter phase is the PLA phase, and the darker phase is the PCL phase).

• The other samples possess a sea-island morphology where the PLA nodules are dispersed in the continuous PCL phase. d→f: The electrical volume conductivity decreases again with the increase of the PLA percentage to reach its minimal value at PLA/10 wt.% GNP. Processing steps:

Double electrical percolation threshold

• Twin-screw extrusion to blend PLA, PCL, and GNP together

• Single-screw extrusion to fabricate calibrated filaments

• Fused filament fabrication to manufacture the samples to be analyzed (compared with the compression molded samples)

• 3D printed sample: Good adhesion between the layers

• 3D printed sample vs compression molded sample: Smaller cocontinuous microstructure existing separately in each deposited layer

Contact information Conclusion References

• This work succeeded for the first time to obtain the double electrical percolation threshold in PLA/PCL/GNP composite systems.

• The influence of several manufacturing processes on the co-continuity of PLA and PCL and consequently the electrical volume conductivity were addressed for the first time.

• Compression molded samples: Wide range of co-continuity that varies between PLA55/PCL45/10 wt.% GNP and PLA70/PCL30/10 wt.% GNP.

Reason: Induced annealing effect

•

  Selective localization of the conductive filler in a single polymer phase • Obtaining the double electrical percolation threshold in poly(lactic acid) (PLA)/ polycaprolactone (PCL)/ Graphene nanoplatelets (GNP) nanocomposites • Studying the influence of different manufacturing processes on the co-continuous microstructure and the electrical conductivity of PLA/PCL/GNP nanocomposites The used processes are: ✓ Compression molding ✓ 3D printing (Fused Filament Fabrication or FFF) preceded by single-screw extrusion In our previous work, a comparison between 3D printing and injection molding influence on the microstructure and therefore the electrical conductivity of PLA/PCL/GNP composites took place [3].

a→d:

  As the percentage of PLA increases from 0 wt.% to 55 wt.% in the composites, the electrical volume conductivity increases due to the greater confinement of the GNP dispersed in the PCL matrix. d: Co-continuous morphologies are observed in the PLA55/PCL45-PLA70/PCL30 range with 10 wt.% GNP. The electrical conductivity is the highest in this zone with the best results for PLA65/PCL35/10 wt.% GNP. A double electrical percolation threshold can be highlighted in this zone.

  3D printed samples: PLA65/PCL35/10 wt.% GNP shows the highest electrical conductivity but still less than that of the compression molded sample Reason: Small co-continuous microstructure present in each deposited layer Corresponding author's information: Nour Alhoda Masarra PhD student at the Ecole des Mines d'Ales (France) [1]: Huang et al., Control of carbon nanotubes at the interface of a co-continuous immiscible polymer blend to fabricate conductive composites with ultralow percolation thresholds, Journal: Carbon, Vol 73, 2014, Pages 267-274 [2]: Chen et al., Balance the electrical properties and mechanical properties of carbon black filled immiscible polymer blends with a double percolation structure, Journal: Composites Science and Technology, Vol 140, 2017, Pages 99-105 [3]: Masarra et al., Fabrication of PLA/PCL/Graphene nanoplatelets (GNP) electrically conductive circuit using the fused filament fabrication (FFF) 3D printing technique,