# Mapping the stability of stellar rotating spheres via linear response theory

S Rozier, J-B Fouvry, P Breen, A Varri, Christophe Pichon, D Heggie

► To cite this version:

S Rozier, J-B Fouvry, P Breen, A Varri, Christophe Pichon, et al.. Mapping the stability of stellar rotating spheres via linear response theory. Monthly Notices of the Royal Astronomical Society, 2019, 487 (1), pp.711-728. 10.1093/mnras/stz1227. hal-02350191

HAL Id: hal-02350191

https://hal.science/hal-02350191

Submitted on 31 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Mapping the stability of stellar rotating spheres via linear response theory

S. Rozier,<fnref n="1" />* J.-B. Fouvry,<fnref n="1" />,<fnref n="2" /> P. G. Breen,<fnref n="3" /> A. L. Varri,<fnref n="4" /> C. Pichon<fnref n="1" />,<fnref n="4" />,<fnref n="5" />
and D. C. Heggie <fnref n="1" />,<fnref n="3" />

Accepted 2019 April 22. Received 2019 April 12; in original form 2019 February 22

## ABSTRACT

Rotation is ubiquitous in the Universe, and recent kinematic surveys have shown that early-type galaxies and globular clusters are no exception. Yet the linear response of spheroidal rotating stellar systems has seldom been studied. This paper takes a step in this direction by considering the behaviour of spherically symmetric systems with differential rotation. Specifically, the stability of several sequences of Plummer spheres is investigated, in which the total angular momentum, as well as the degree and flavour of anisotropy in the velocity space are varied. To that end, the response matrix method is customized to spherical rotating equilibria. The shapes, pattern speeds and growth rates of the systems' unstable modes are computed. Detailed comparisons to appropriate N-body measurements are also presented. The marginal stability boundary is charted in the parameter space of velocity anisotropy and rotation rate. When rotation is introduced, two sequences of growing modes are identified corresponding to radially and tangentially biased anisotropic spheres, respectively. For radially anisotropic spheres, growing modes occur on two intersecting surfaces (in the parameter space of anisotropy and rotation), which correspond to fast and slow modes, depending on the net rotation rate. Generalized, approximate stability criteria are finally presented.

**Key words:** diffusion – gravitation – galaxies: kinematics and dynamics – galaxies: nuclei.

## 1 INTRODUCTION

Most stellar systems are rotating. While this has been clear for spiral galaxies since the beginning of extragalactic spectroscopy, recent kinematic surveys have shown that early-type galaxies and globular clusters can also be partially centrifugally supported.

Globular clusters have offered an essential empirical probe to develop and test theories of stellar dynamics. The long-term evolution of such self-gravitating systems generically comprises two stages. First, on a few dynamical times, as a result of strong potential fluctuations, the cluster can undergo a violent (collisionless) relaxation (Lynden-Bell 1967), allowing it to reach a steady state. Then, on longer time-scales, the cluster will slowly explore subsequent thermodynamical equilibria, as a result of two-body relaxation sourced by finite-N fluctuations (Binney & Tremaine 2008). During the first stage of evolution, the globular cluster is insensitive to its finite number of constituents, and can therefore be treated as a collisionless system. It is also during this stage that the system can develop linear instabilities, that deeply rearrange the

system's orbital structure, should the violent relaxation stage have left it in an unlikely (low entropy) configuration. This is where lies the interest of this paper.

The analysis of the linear stability of spherically symmetric stellar systems has already been the subject of numerous investigations. For example, the Doremus-Feix-Baumann theorem (Binney & Tremaine 2008) states that radial modes of an ergodic spherical model are all stable if the gradient of the system's distribution function (DF) with respect to the energy is negative, while systems with a positive gradient with respect to the angular momentum tend to be unstable (Tremaine 2005). Similarly, systems that support too many radial orbits tend to undergo the so-called radial orbit instability (ROI, see e.g. Polyachenko & Shukhman 1981; Saha 1991; Weinberg 1991). However, all these investigations were limited to configurations having a spherically symmetric velocity distribution (i.e. a DF which is invariant under rotation in configuration space). This study sets out to lift this restriction and investigate the linear stability of rotating, anisotropic equilibria, i.e. systems with a non-zero total angular momentum.

This theoretical endeavour is particularly timely, for a number of reasons. The interest in the role of angular momentum in the

<footnote n="1">¹Institut d'Astrophysique de Paris, UPMC Univ. Paris 06, (UMR7095), 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="2">²Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540, USA</footnote>

<footnote n="3">³School of Mathematics and Maxwell Institute for Mathematical Sciences, University of Edinburgh, Kings Buildings, Edinburgh EH9 3FD, UK</footnote>

<footnote n="4">⁴Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="5">⁵Korea Institute for Advanced Study (KIAS), 85 Hoegiro, Dongdaemun-gu, Seoul, 02455, Republic of Korea</footnote>

<footnote>* E-mail: rozier@iap.fr</footnote>
evolution of globular clusters is now peaking, as a result of a variety of new empirical evidence. This includes recent astrometric measurements by Hubble Space Telescope (HST, Bellini et al. 2017), Gaia DR2 (Bianchini et al. 2018; Sollima, Baumgardt & Hilker 2019), as well as spectroscopic surveys, such as the Multi-Instrument Kinematic Survey (Ferraro et al. 2018) and dedicated programmes with the Multi Unit Spectroscopic Explorer (Kamann et al. 2018). Similarly, line-of-sight measurements of nuclear star clusters are also revealing them as fast rotators (Seth et al. 2008; Feldmeier-Krause et al. 2017). Finally, on larger scales, the classical V/σ versus ellipticity diagram (Davies et al. 1983) is now rich in spheroidal, rotation-supported early-type galaxies (see e.g. Davies 2011 for a summary), following major surveys such as SAURON (Bacon et al. 2001) and ATLAS3D (Cappellari et al. 2011). While these systems are often modelled with Jeans and Schwarzschild techniques, it is important to explore the stability properties of the equilibria identified in these studies.

The main goal of this work is to eventually extend to stellar systems the classical knowledge of the stability of uniformly rotating fluid dynamical systems (for a summary of the theory of the ellipsoidal figures of equilibrium, see Chandrasekhar 1969), and its subsequent developments into the realm of differentially rotating polytropes, which have been found to exhibit instability at surprisingly low values of the total angular momentum (e.g. see Pickett, Durisen & Davis 1996; New & Shapiro 2001; Centrella et al. 2001; Shibata, Karino & Eriguchi 2002; Ou & Tohline 2006 and, more recently, Galeazzi, Yoshida & Eriguchi 2012; suggestions for a possible dynamical interpretation in which corotation points play a primary role have been offered, among others, by Watts, Andersson & Jones 2005; Passamonti & Andersson 2015; Yoshida & Saijo 2017).

The present study is therefore part of a concerted effort to systematically explore the theoretical dimension of the newly appreciated 'kinematic richness' (i.e. any deviations from the traditional assumptions of isotropy in the velocity space and absence of internal rotation), which is empirically emerging in a variety of spheroidal stellar systems. Such a richness has implications on (i) the construction of appropriate self-consistent equilibria (e.g. see the fresh efforts by Binney 2014; Evans et al. 2015; Sanders & Evans 2015; Posti et al. 2015 and, for a summary of recent developments with applications to globular clusters, Varri 2016, and references therein), (ii) their dynamical and secular (collisionless) behaviour, as well as (iii) the long-term (collisional) evolution of stellar systems (for an analysis of the effects of 'primordial' velocity anisotropy and differential rotation on spherically symmetric collisional systems, see Breen, Varri & Heggie 2017; Breen et al. 2019a, respectively).

With the exception of some pioneering studies on the counterparts of uniformly rotating polytropes (see especially Vandervoort 1980), very few investigations of the stability properties of spheroidal, differentially rotating stellar systems have been conducted so far, and mostly by means of numerical approaches (see e.g. Barnes, Goodman & Hut 1986; Kuijken & Dubinski 1994; van der Marel, Sigurdsson & Hernquist 1997; Varri 2012), partly with the hope of exploring the validity of the classic (Ostriker & Peebles 1973) stability criterion in that context. It is therefore opportune to explore the role of angular momentum in the early evolution of spheroidal stellar systems and investigate the possible interplay between velocity anisotropies and internal rotation.

This paper is organized as follows. Section 2 presents the equilibrium model which we consider as initial state, in particular its DF with velocity anisotropies, and how internal rotation is introduced

(through Lynden-Bell's demon). With the aim of characterizing the linear stability of these systems, Section 3 presents an extension of the matrix method for rotating spherical systems, applies this theory to rotating, anisotropic Plummer spheres, and provides also a direct comparison with N-body simulations. Section 4 extends this analysis to the whole 2D parameter space of velocity anisotropy and rotation, and investigates their respective roles in sourcing dynamical instabilities. We finally discuss our results in Section 5, and conclude in Section 6.

## 2 DESCRIBING ROTATING SPHERES

Let us start by describing in detail the equilibrium model adopted for the initial state of the system, focusing on its two-parameter DF that encodes both the velocity anisotropies and the system's total rotation.

### 2.1 Angle-action coordinates

The mean stellar system is always assumed to be spherically symmetric, and characterized by a mean radial potential $\psi(r)$, with $(r, \theta, \phi)$ the spherical coordinates. Throughout this study, we take the mean potential to be a Plummer potential (Plummer 1911),

$$ \psi(r) = -\frac{GM_{\text{tot}}}{R_s} \frac{1}{\sqrt{1 + (r/R_s)^2}}, \quad (1) $$

and work in the dimensionless units $G = R_s = M_{\text{tot}} = 1$, with $M_{\text{tot}}$ the system's total mass. A spherically symmetric potential guarantees the existence of angle-action coordinates, $(\theta, \mathbf{J})$, (Binney & Tremaine 2008), with

$$ \mathbf{J} = (J_1, J_2, J_3) = (J_r, L, L_z), \quad (2) $$

standing respectively for the radial action, the norm of the angular momentum, and its projection along the vertical axis. The orbit's inclination is naturally defined as $\cos(i) = L_z/L$, with $0 \le i \le \pi$. According to the Jeans theorem (Jeans 1915), the mean-field distribution of the system is fully characterized by its DF, $F(\mathbf{J})$, normalized as $\int d\theta d\mathbf{J} F(\mathbf{J}) = M_{\text{tot}}$. Because the system is spherically symmetric, the orbital frequencies, $\Omega(\mathbf{J}) = \partial H/\partial \mathbf{J}$, satisfy the degeneracy condition $\Omega_3(\mathbf{J}) = 0$, i.e. stars orbit in a fixed orbital plane. Similarly to Hamilton et al. (2018), this already entices us to define the reduced actions $\tilde{\mathbf{J}} = (J_r, L)$, and the associated reduced frequencies $\tilde{\Omega} = \partial H(\tilde{\mathbf{J}})/\partial \tilde{\mathbf{J}}$. Note that spherical systems with a spherically symmetric velocity distribution satisfy $F = F(\tilde{\mathbf{J}})$.

### 2.2 Introducing velocity anisotropies

As in the investigation conducted by Breen et al. (2017), in order to study possible effects associated with velocity anisotropies, we consider the family of DFs for the Plummer sphere proposed by Dejonghe (1987)

$$ F_0(q, E, L) = \frac{3\Gamma(6-q)}{2(2\pi)^{5/2}} (-E)^{\frac{7-q}{2}} \\ \times \mathcal{F} \left( 0, q/2, \frac{9}{2} - q, 1; -\frac{L^2}{2E} \right), \quad (3) $$

with $q$ controlling the flavour and degree of velocity anisotropy of the system, and $\Gamma$ the gamma function. Here, $\mathcal{F}$ is expressed in
<figure><img src="image_3.png" /><figcaption>Figure 1. Illustration of the radial dependence of $\langle v_{\phi}(r) \rangle$ (top panel) and $\beta_{\phi}(r)$ (bottom panel) for various velocity anisotropies (via $q$, indicated by the shade of colour) and total angular momentum (via $\alpha$, indicated by the line style). Here, $\alpha$ controls $\langle v_{\phi}(r) \rangle$, while $q$ controls $\beta_{\phi}$, noticeably in the outer region of the cluster. We note that for $(\alpha, q) = (1, -6)$ a significant fraction of the radial support is centrifugal, whereas $\sigma_{\phi}$ remains comparatively small. The sphere has therefore a large reservoir of kinetic energy, which can be converted into thermal energy through an instability.</figcaption></figure>

terms<fnref n="1" /> of the hypergeometric function $_2F_1$

$F(a, b, c, d; x)$

$$ = \begin{cases} \frac{x^a {}_2F_1(a+b, 1+a-c; a+d; x)}{\Gamma(c-a)\Gamma(a+d)}, & \text{if } x \le 1; \\ \frac{1}{{}_x^b {}_2F_1(a+b, 1+b-d; b+c; \frac{1}{x})} \frac{r^2}{\Gamma(d-b)\Gamma(b+c)}, & \text{if } x \ge 1. \end{cases} \quad (4) $$

As illustrated in Fig. 1, the parameter $q$ controls the flavour of the velocity anisotropy in the system, with $q < 0$ (respectively, $q > 0$) associated with tangential (respectively, radial) velocity anisotropies. One of the advantages of this DF is the simple forms taken by its velocity dispersions in all directions. Indeed, one has

$$ \sigma_r^2(r) = \frac{(1 + r^2)^{-1/2}}{6 - q}, $$

$$ \sigma_{\theta}^{2}(r) = \sigma_{\phi}^{2}(r) = \sigma_{r}^{2}(r) \left( 1 - \frac{q}{2} \frac{r^2}{1 + r^2} \right), \quad (5) $$

leading to the anisotropy parameter:

$$ \beta_{\theta}(r) = 1 - \frac{\sigma_{\theta}^{2}(r)}{\sigma_{r}^{2}(r)} = \frac{q}{2} \frac{r^2}{1 + r^2}. \quad (6) $$

We note that after introducing rotation as below, this relation still holds, as long as one uses equation (6) to define $\beta_{\theta}$, ignoring the dispersion in $\phi$.

## 2.3 Introducing rotation

Following Lynden-Bell (1960), a simple way of introducing rotation in a spherically symmetric system is to modify its DF as<fnref n="2" />

$$ F(\alpha, q, \mathbf{J}) = F_0(q, \tilde{\mathbf{J}}) + \alpha F_0(q, \tilde{\mathbf{J}}) \text{Sign}(L_z), \quad (7) $$

with the sign function $\text{Sign}(x)$ equal to $-1$ for $x < 0$, then $0$ for $x = 0$, and $1$ for $x > 0$. In equation (7), $|\alpha| \le 1$ is a dimensionless parameter that controls the amount of rotation in the system. We note that $F(\alpha = 0, q, \mathbf{J}) = F_0(q, \tilde{\mathbf{J}})$, corresponding to the case without rotation, and that the normalization of $F$ is guaranteed by that of $F_0$, since $\int_{-L}^{L} dL_z \text{Sign}(L_z) = 0$. The DF of equation (7) does not change the distribution of the orbits' shapes, but rather the distribution of their orientations. For each orbit (i.e. each $\tilde{\mathbf{J}}$), the probability to have $L_z > 0$ is multiplied by $1 + \alpha$, compared to the non-rotating DF. In the limiting case $\alpha = 1$, all particles have $L_z > 0$.

Since a more extensive characterization of the properties of these configurations will be provided in Breen et al. (2019a), here we just wish to highlight a few selected features.<fnref n="3" /> In Fig. 1, we illustrate how $\alpha$ creates some mean azimuthal motion in an otherwise non-rotating system. We also note that $\alpha$ induces an anisotropy in the azimuthal motions, by breaking the symmetry in the two tangential directions, as given by the azimuthal anisotropy parameter

$$ \beta_{\phi}(r) = 1 - \frac{\sigma_{\phi, \alpha}^{2}(r)}{\sigma_{r}^{2}(r)}. \quad (8) $$

The bottom panel of Fig. 1 also shows how $\alpha$ biases the azimuthal anisotropy towards radial motions by increasing the mean azimuthal velocity. Similarly, increasing $q$ increases the mean azimuthal velocity in the central region of the cluster and decreases it in the outer region. For $q \ne 0$, a radial variation of $\beta_{\phi}$ is also introduced in the outer region of the cluster. Moreover, increasing $q$ decreases $\beta_{\phi}$. We note that by construction $\langle v_{\phi}(r) \rangle \propto \alpha$. Finally, we note that even for $q = -6$, when $\alpha = 1$, $\sigma_r > \sigma_{\phi}$.

In Fig. 2, we illustrate the probability distribution function (PDF) of $(v_r, v_{\phi})$ for various values of $\alpha$ and $q$. We note that the ansatz given by equation (7) introduces a sharp discontinuity in the distribution of stars in the $v_{\phi}$-direction as soon as $\alpha \ne 0$. It is possible that this discontinuity plays a role in driving the instabilities, as discussed further in Section 5.2.

# 3 CHARACTERIZING UNSTABLE MODES

## 3.1 Matrix method and rotating spheres

The stability properties of anisotropic, spherical equilibria have already been the object of many studies (Polyachenko & Shukhman

<footnote n="1">¹This formulation slightly simplifies Dejonghe's by suppressing vanishing $\Gamma(q/2)$ terms, which are indeterminate when $q$ is an even, negative integer.</footnote>

<footnote n="2">²We note that any odd function of $L_z$ could be used to generalize our choice of Sign($L_z$) to introduce rotation.</footnote>

<footnote n="3">³An implementation of the construction of such anisotropic, differentially rotating spherical models is already freely available at https://github.com/pgbreen/PlummerPlus.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Number density of stars in the ($v_r$, $v_\phi$)-space at the constant radius $r=1$. Whenever $\alpha \neq 0$, the PDF is discontinuous on the line $v_\phi = 0$. As $q$ increases, the spread in $v_\phi$ decreases, and the system becomes more radially biased. As $\alpha$ increases, stars are swapped towards positive $v_\phi$.</figcaption></figure>

1981; Saha 1991; Weinberg 1991; Polyachenko & Shukhman 2015). Below, we follow the notations from Hamilton et al. (2018).

The linear stability of a long-range interacting integrable system can generically be characterized by its response matrix $\hat{M}(\omega)$ (Binyey & Tremaine 2008), which in the case of 3D spherical systems is given by

$$ \hat{M}_{pq}(\omega) = (2\pi)^3 \sum_{n \in \mathbb{Z}^3} \int d\mathbf{J} \frac{\mathbf{n} \cdot \partial F / \partial \mathbf{J}}{\omega - \mathbf{n} \cdot \mathbf{\Omega}(\mathbf{J})} [\psi_n^{(p)}(\mathbf{J})]^* \psi_n^{(q)}(\mathbf{J}). \quad (9) $$

A system is said to be linearly unstable if there exists a frequency $\omega = \omega_0 + i\eta$ (with $\eta > 0$), such that $\hat{M}(\omega)$ admits an eigenvalue equal to 1. The instability criterion for the sphere therefore reads

$$ \det[\mathbf{I} - \hat{M}(\omega_0 + i\eta)] = 0. \quad (10) $$

In that case, the system supports an unstable mode of oscillation frequency $\omega_0$ and growth rate $\eta$. To capture the system's self-gravitating response, equation (9) involves a biorthogonal basis of potentials and densities, $\psi^{(p)}(\mathbf{x})$ and $\rho^{(p)}(\mathbf{x})$, whose conventions are given explicitly in Appendix A. The same appendix details how the response matrix associated with the rotating DF from equation (7) is computed, and shows that

$$ \hat{M}(\alpha, \omega) = \hat{M}_0(\omega) + \alpha \hat{M}_1(\omega), \quad (11) $$

where $\hat{M}_0(\omega)$ corresponds to the system's response in the absence of rotation, while $\alpha \hat{M}_1(\omega)$ captures the effects associated with rotation. Equation (11) is the master equation characterizing its linear stability.<fnref n="4" /> Because of its reduced number of symmetries, we

also note that the matrix $\hat{M}_1(\omega)$ is less sparse than the zero-rotation matrix $\hat{M}_0(\omega)$, and hence the numerical task is more cumbersome. In Appendix B, we detail our numerical implementation of the computation of $\hat{M}(\alpha, \omega)$.

Before applying the previous machinery to the study of the linear stability of a wide range of equilibria [i.e. different $(\alpha, q)$], we will start by considering one particular case. Following equation (7), we focus on the fiducial model $(\alpha, q) = (1, -6)$, i.e. a model with the maximum rotation and with tangential velocity anisotropies, as already illustrated in Fig. 1. From now on, our search is restricted to unstable modes with the pattern number $m=2$, i.e. bars and two-armed spirals. This is motivated by previous results from N-body simulations, which were conducted in continuation of Breen et al. (2017) and will appear in a separate study (Breen et al. 2019b), which suggested that bisymmetric modes dominate.

## 3.2 Matrix method methodology

In order to investigate the stability properties of the fiducial model, the numerical methods detailed in Appendices A and B are implemented. Following Pichon & Cannon (1997), unstable modes are identified using Nyquist contours. To do so, for a given value of $\eta > 0$, one constructs the curve $\omega_0 \mapsto \det[\mathbf{I} - \hat{M}(\omega_0 + i\eta)]$. If ever this contour encloses the origin, then the cluster supports an unstable mode with a growth rate larger than $\eta$. Figure 3 illustrates the Nyquist contours associated with our fiducial model.

The set of numerical parameters of this computation is detailed in Appendix C1, where we also explore the convergence of the matrix method with respect to these various numerical parameters. In order to find the instability with the largest growth rate, we draw Nyquist contours for progressively decreasing values of $\eta$, from a value at which the Nyquist contour does not enclose the origin, down to the first crossing point. For the fiducial model, we identify an unstable mode with frequency $(\omega_0, \eta) \simeq (0.547, 0.047)$. We note that as one keeps decreasing the value of $\eta$, additional subdominant modes may be found. This will be of particular interest in Section 4.3, when considering radially biased systems.

Once the unstable mode has been identified, the response matrix's eigenvector associated with the eigenvalue nearly equal to 1 can be determined, which allows us to obtain the mode's shape via

$$ \rho(r, \vartheta, \phi) = \mathrm{Re} \left( \sum_{\rho} X_{\rho} \rho^{(\rho)}(r, \vartheta, \phi) \right), \quad (12) $$

with $X$ the eigenvector. In Fig. 4, we represent in the equatorial plane the density of the mode of the fiducial model. In Fig. 5, we illustrate the 3D density isocontours of the same unstable mode. We note that the density response exhibits a large-scale spiral pattern extending beyond corotation, and that the mode's thickness in the vertical direction is significant.

## 3.3 N-body methodology

In order to validate the linear stability analysis further, N-body simulations were carried out for the same fiducial model, with $N=2^{17}$ particles, integrated for a time interval of 100 Hénon Units (Hénon 1971). The simulations were performed using GYRFALCON (Dehnen 2002). To characterize the frequencies of the unstable modes, we Fourier-analysed the space density as follows. For a given pattern number $m=2$, we define $A_m(t) = C_m(t)/C_0(t)$, where

$$ C_m(t) = \int_0^{+\infty} dr \, r \int_{-\infty}^{+\infty} dz \int_0^{2\pi} \frac{d\phi}{2\pi} \rho(r, z, \phi, t) e^{-im\phi}, \quad (13) $$

<footnote n="4">⁴We note that in the slow rotating rate regime, i.e. $\alpha \ll 1$, seeking eigenvalues for equation (11) is formally identical to the corresponding problem addressed by perturbation theory in quantum mechanics. In particular, given the knowledge of the full spectrum of $\hat{M}_0(\omega)$, the shift in eigenvalues and eigenvectors associated with the introduction of $\alpha \ll 1$ can be computed perturbatively.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Nyquist diagram of the fastest growing mode of the model $(\alpha, q) = (1, -6)$, as given by the curves $\omega_0 \mapsto \det[\mathbf{I} - \hat{\mathbf{M}}(\omega_0 + i\eta)]$. The top panel represents the Nyquist diagram without renormalization. Following Pichon & Cannon (1997), in the bottom panel, a logarithmic scaling is applied, $re^{i\theta} \mapsto 1/3\log(1 + 10^3r)e^{i\theta}$, to highlight the structure of the contour around the origin and make sure that the first growing mode is indeed selected. Different colours are associated with different growth rates, while the changes in $\omega_0$ along each curve are indicated with the grey dashed lines. This diagram shows that the cluster supports an unstable mode with a growth rate, $0.045 \le \eta \le 0.049$, and oscillation frequency, $0.547 \le \omega_0 \le 0.550$.</figcaption></figure>

with $\rho(R, z, \phi, t)$ the simulation's instantaneous density (i.e. a sum of N delta functions). Should the cluster support an unstable mode, one can in principle extract its frequencies from the numerically measured function $t \mapsto A_m(t)$, as one has

$$ \frac{d \log(|A_m(t)|)}{dt} = \eta; \quad \frac{d \operatorname{Arg}(A_m(t))}{dt} = \omega_0. \qquad (14) $$

In practice, measuring the growth rate of an instability in N-body data has two major difficulties. The first is noise due to the discrete nature of the system. The second is that the N-body system cannot sustain an exponential growth indefinitely, and the mode must saturate at some point. The transition from exponential growth to saturation is not easy to model, and may have some dependence on the phase space DF. The result is that this transition can possibly result in a bias such that growth rates based on an exponential fit may systematically underestimate the true growth actually experienced by the mode in its early linear phase.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Radial dependence in the equatorial plane of the density of the unstable mode recovered in Figs. 3 and 6. The top panel is the linear theory's prediction, and the bottom panel the measurement from N-body simulations. The black dashed circle represents the scale radius of the background Plummer potential, and the red circle the mode's corotation radius. Isocontours are logarithmically sampled, with a factor ~1.6 between two consecutive lines, so that the outermost line has a density $1.5 \times 10^4$ times smaller than the peak value. The similarity between the two results is striking.</figcaption></figure>

In order to automatically detect and measure growth rates in a grid of N-body simulations, we use an iterative scheme which involves fitting an exponential curve to the amplitude $|A_m(t)|$ of the mode, first over the full duration of the simulation, but then reducing the time interval of the fit, until the maximum absolute difference between the fit and data (the normalized amplitude) falls below a pre-determined value of tolerance. The choice of tolerance
<figure><img src="image_7.png" /><figcaption>**Figure 5.** 3D isosurfaces of the density of the unstable mode recovered in Fig. 3, and their intersection with the equatorial and meridional planes.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Illustration of the time-dependence of $|A_m(t)|$ (left-hand panel) and Arg($A_m(t)$) (right-hand panel) for the same cluster model as in Fig. 3. Following equation (14), linear fits are also represented. We recovered ($\omega_0$, $\eta$) = (0.57, 0.022) that should be compared with the measurements of Fig. 3.</figcaption></figure>

should be small enough such that the time interval is reduced when there is a poor fit to a mode (because of its saturation) but large enough so that the interval is not reduced merely by the presence of noise. Based on a grid of 300 different realizations of the fiducial model ($N = 2^{17}$, $q = -6$, and $\alpha = 1$), a value of 0.009 was found to produce robust measurements of the growth rate $\eta$. Similarly, $\omega_0$ is determined by a linear fit to Arg($A_m(t)$) (mod $2\pi$ to account for complete revolutions) over the same time interval. Figure 6 illustrates such measurements of $\eta$ and $\omega_0$ for the fiducial model.

Having characterized the frequencies of the unstable modes, one can use the same numerical simulations to extract the mode's shape. To do so, the simulation's particles are projected on to the same $\rho^{(p)}(r, \varphi)$ basis as for the matrix method, to reconstruct the mode using equation (12). This is illustrated in the bottom panel of Fig. 4.

In Appendix C2, we briefly check the convergence of the N-body runs with respect to their various control parameters.

The comparison between linear theory and N-body simulations in the fiducial model shows a ~15 per cent discrepancy in the measured value of $\eta$, and a ~3 per cent discrepancy in $\omega_0$. Moreover, as highlighted in Fig. 4, the modes' shapes are in close agreement. This allows us to quantitatively validate our implementation of the matrix method as well as our extension of the linear stability analysis to rotating systems.

# 4 CHARTING THE INSTABILITIES

Having found good agreement between the matrix method and the direct N-body approach, we now set out to investigate the properties of the unstable modes in the whole $(\alpha, q)$-space, to account for the contributions from both velocity anisotropy and total angular momentum. As in the previous sections, we will focus on the three characteristics of unstable modes, namely the growth rate $\eta$, the oscillation frequency $\omega_0$, and the mode's shape.

## 4.1 Instability in the $(\alpha, q)$-space

We reproduced the methodology of Section 3.2 to chart the $(\alpha, q)$-space using the linear matrix method, as well as N-body simulations. The resulting frequencies are illustrated in Fig. 7.

In that figure, one can straightforwardly identify two types of instabilities. The first one, corresponding to the pyramid on the top of the $\eta(\alpha, q)$ panel ($q < 0$, $\alpha \gtrsim 1/2$) appears to be specific to tangentially biased and rotating systems. These instabilities are studied in detail in Section 4.2. A second type of instability can also be noted in the bottom right part of the $\eta(\alpha, q)$ panel ($q > 0$). These instabilities appear to be specific to radially biased systems. We investigate these instabilities in Section 4.3. The two right-hand panels of Fig. 7 show the N-body counterparts of the measurements using the matrix method. Again, the quantitative matchings between the left- and right-hand panels validate our generalization of the matrix method. On all these panels, for both $\eta$ and $\omega_0$, the lower blue plane corresponds to a growth rate threshold, $\eta = 0.01$, below which we stopped the search for unstable modes, deeming that the accuracy of the methods is not sufficient anymore, at least with the choice of accuracy parameters made in this work (see Appendix C).<fnref n="5" />

Having validated the methodology against N-body simulations, we may now carry out a more thorough exploration of the $(\alpha, q)$-space using only the matrix method, successively for tangentially ($q < 0$) and radially biased systems ($q > 0$).

## 4.2 Instability of tangentially biased systems: $q < 0$

In Fig. 8, we illustrate cuts of constant $q$ in the $\eta(\alpha, q)$ and $\omega(\alpha, q)$ surfaces in the $q < 0$ region. In the relevant region, the instability surface $(\alpha, q) \to \omega_0$, $\eta$ can be fitted globally by the simple polynomial expressions

$$ \begin{aligned} \omega_0 &= \frac{11\alpha}{15} - \frac{q}{66} - \frac{2}{9}, \\ \eta &= \frac{29\alpha^2}{50} - \frac{q^2}{8333} - \frac{2q\alpha}{111} + \frac{q}{103} - \frac{47\alpha}{50} + \frac{37}{100}. \end{aligned} \quad (15) $$

<footnote n="5">⁵Indeed, we found that the matrix method, while potentially optimal to identify linear growth rates, becomes cumbersome in the small growth rate regime (see also Merritt 1999).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Illustration of the dependence of the growth rates $\eta$ (top panels) and oscillation frequencies $\omega_0$ (bottom panels), as a function of the cluster's parameters $(\alpha, q)$, measured using the response matrix method (left-hand panels) and $N$-body methods (right-hand panels). We searched for unstable modes on a $(\alpha, q)$-grid composed of the locations $\alpha = 0, 0.2, 0.4, 0.6, 0.8$, and $1$ and $q = -16, -12, -6, -2, 0, 1$, and $2$. The blue plane represents our growth rate threshold, $\eta = 0.01$, for mode selection. At this level of accuracy, it appears that rotation introduces a new island of instability in the top left region of each panel (corresponding to $q < 0, \alpha > 1/2$) and an extension of the known ROI instability island on the right (near $q \sim 2$). The comparison of the matrix method's and $N$-body measurements highlights again the reliability of our implementation of the linear theory.</figcaption></figure>

These fitting formulae, as well as the fitting formulae we present for other instability surfaces, were obtained by applying a least-squares algorithm to different low-order polynomial forms, and choosing the fit that both guarantees a reasonable accuracy with respect to the data and a simple functional form. The cuts of the present fitting formulae are represented in Fig. 8, showing their level of accuracy.

In Fig. 9, we represent the $\eta(\alpha, q)$ surface for tangentially biased systems. We emphasize that in this regime, both $\eta$ and $\omega_0$ increase with $\alpha$, and increase as $q$ decreases. We also note that, so far, these instabilities appear to emerge only in the presence of non-vanishing angular momentum (i.e. they are not present in configuration with $\alpha \to 0$). We also emphasize that their specific behaviour depend on the adopted choice for the angular momentum distribution (i.e. as resulting from the application of Lynden-Bell's demon), as discussed in Section 5.2. In Fig. 9, we also show the shape of the

modes at different locations of the instability surface. The absence of any significant variation in the modes' shape suggests that the underlying instability mechanism remains the same.

### 4.3 Instability of radially biased systems: $q > 0$

In Fig. 10, we illustrate cuts of constant $q$ in the $\eta(\alpha, q)$ and $\omega(\alpha, q)$ surfaces in the $q > 0$ region. Similarly, in Fig. 11, we explore the $\eta(\alpha, q)$ surface, and represent the shapes of the associated unstable modes. The shapes of the modes suggest in this regime the existence of two types of instability, depending on the region considered.

The first type of instability dominates for slow-rotating systems ($\alpha \lesssim 1/2$). It is represented with the solid lines in Fig. 10, and the dark grey surface in Fig. 11. It can be identified as the continuation in the rotating regime of the standard ROI (Polyachenko & Shukhman 1981). This type will be referred to hereafter as 'Fast-ROI', since the
<figure><img src="image_10.png" /><figcaption>Figure 8. Illustration of the growth rate (top panel) and oscillation frequency (bottom panel) as a function of $\alpha$ for fixed values of $q$, in the $q < 0$ region. The grey region represents the threshold for the measurement of instabilities. At this level of accuracy, the $\eta$ and $\omega_0$ surfaces can be fitted globally by quadratic functions of $\alpha$ and $q$, represented with the thin lines.</figcaption></figure>

oscillation frequency of the instability is large. For that particular instability, the modes evolve from a bar-like shape around $q=0$ to a loosely wound trailing spiral around $q=2$, as shown in Fig. 11.

The second type of instability dominates for fast-rotating systems ($\alpha \gtrsim 1/2$). It is represented with the dashed lines in Fig. 10, and the light grey surface in Fig. 11. We note that the oscillation frequencies of these modes are small, and we therefore refer to this instability as ‘Slow-ROI’. As illustrated in Fig. 11, for that second instability regime, the modes present a barred shape, with a more wound trailing spiral departing from the bar, that qualitatively differs from the shape of the Fast-ROI modes.

Both instability surfaces ($\alpha, q) \to \omega_0, \eta$ can be fitted globally by simple polynomial functions as

$$ 
\begin{aligned}
\omega_0 &= \frac{4\alpha}{11} - \frac{q}{34} + \frac{1}{38}, & \eta &= \frac{q}{76} + \frac{11\alpha}{100} - \frac{q\alpha}{109} - \frac{7\alpha^2}{50} - \frac{1}{66}, \\
\omega_0 &= \frac{2\alpha}{25} - \frac{q}{123} + \frac{1}{39}, & \eta &= \frac{q}{52} + \frac{\alpha}{23} - \frac{1}{28},
\end{aligned} 
\quad (16) 
$$

for the Fast- and the Slow-ROI, respectively.

In Fig. 12, we finally illustrate isocontours of the growth rates of the dominant instabilities in the $(\alpha, q > 0)$-space. This figure clearly highlights the transition between the two regions of instability, Fast-ROI versus Slow-ROI.

## 5 DISCUSSION

The previous sections identified various islands of instability in the two parameters $(\alpha, q)$ for the model from Section 2. Such a

charting is costly for rotating spheres, because the broken symmetry makes the matrix problem larger, as azimuthal harmonics couple via gravity (see Appendix A). Moreover, the previous charting is clearly particular to our choice of potential, as well as our parametrization of the DF and the angular momentum distribution. However, it nonetheless displays some diversity, which hopefully should be of general character. Let us therefore remap the stability chart in terms of physically motivated properties of rotating spheres, before commenting on the origin and the topography of the instability islands.

### 5.1 Remapping the instability chart

In order to rephrase our results in a more general framework, we shall define macroscopic kinematic parameters of the sphere via the mass-weighted radial average of $\beta_\theta$ and $L_\alpha$ respectively, using equations (6) and (7). They read

$$ 
\begin{aligned}
\langle \beta_{\theta} \rangle &\equiv \frac{1}{M_{\text{tot}}} \int dr \, 4\pi r^2 \rho(r) \beta_{\theta}(r) = \frac{3}{10} q, \\
\langle L_z \rangle &= \alpha \sqrt{G M_{\text{tot}} R_s} l(q), \quad \text{with } l(q) \approx \frac{3}{5} - \frac{10}{9} \left(\frac{69}{14} - q\right)^{-1}, (17)
\end{aligned}
$$

where $l(q)$ is a good fit decreasing monotonically from 0.6 to 0.22 as $q$ increases. From equation (17), we note that to zeroth order, $q$ and $\alpha$ are direct proxy for $\langle \beta_\theta \rangle$ and $\langle L_z \rangle$. Using this reparametrization, we present in Fig. 13 a remapping of the stability manifolds in the $(\langle \beta_\theta \rangle, \langle L_z \rangle)$-space.

As expected the qualitative differences between Figs 7 and 13 are minor. The following simple fits represent the three manifolds

$$ 
\begin{aligned}
\omega_0 &= \frac{10L}{7} - \frac{\beta}{43} - \frac{1}{7}, & \eta &= \frac{117L^2}{55} - \frac{\beta^2}{1623} - \frac{L\beta}{26} - \frac{141L}{89} \\
& & + \frac{\beta}{107} + \frac{14}{47}, \\
\omega_0 &= L, & \eta &= \frac{11L}{21} - \frac{63L^2}{38} + \frac{2\beta}{17} - \frac{18L\beta}{37} - \frac{\beta^2}{20} - \frac{1}{28}, \\
\omega_0 &= \frac{L}{3}, & \eta &= \frac{L}{7} + \frac{\beta}{10} - \frac{1}{21},
\end{aligned} 
\quad (18) 
$$

respectively for negative $q$, Fast-ROI, Slow-ROI, using the shortening notations $\beta = \langle \beta_\theta \rangle$ and $L = \langle L_z \rangle$. From equation (18), the Ostriker-Peebles criterion (Ostriker & Peebles 1973) can be generalized by requesting that $\eta$ is above a given threshold $\eta_{\text{crit}}$. For instance, for the Slow-ROI instability, it reads $\langle L_z \rangle/7 + \langle \beta_\theta \rangle/10 \ge \eta_{\text{crit}} + 1/21$.

### 5.2 The role of the DF's discontinuity

It is generally now accepted that DFs produced by violent or secular relaxation need not be very smooth in action space (e.g. Fouvry et al. 2015). Yet, typically, these discontinuities are not of the form determined by the application of Lynden-Bell's demon from equation (7) which may partially impact artificially the distribution of orbits near $L_\alpha = 0$, as can be seen in Fig. 2. This may prevent other instabilities from playing a role within the top left part of the chart from Fig. 7, though it should be noted that clusters with $q \sim 0$ and $\alpha \sim 1$ are not unstable, at least according to the criterion adopted here.

In order to check if the specific discontinuity in equation (7) is responsible for the existence of the modes, we have applied the
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Illustration of the various mode shapes as one explores the plane of instability in the $q < 0$ region. The black dashed circle represents the scale radius of the Plummer potential. The grey arrow indicates the direction of rotation of both the global rotation and the instability patterns. In addition to matching the colour of the mode with the points in the central figure, the value of $(\alpha, q)$ is written inside each circle. With increasing $\alpha$ and decreasing $q$, the spiral patterns get more and more wound.</figcaption></figure>

matrix method to a smoothed modification of the DF, in which the Sign function is replaced by the expression

$$ \frac{\operatorname{erf}(a L_z / L)}{\operatorname{erf}(a)}, $$

where $a$ is a positive real parameter. In the limit $a \to 0$, this tends to the function $L_z/L$, whereas in the limit $a \to \infty$, it tends pointwise to $\operatorname{Sign}(L_z)$. For the standard case $(\alpha, q) = (1, -6)$, the growth rate $\eta$ increases monotonically from 0.0060 when $a = 0.1$ to 0.045 when $a = 100$. With the discontinuous DF, the corresponding result is $\eta = 0.045$ (Table C1), in good agreement with the result for $a = 100$. The pattern speed $\omega$ increases from 0.30 to 0.54 as $a$ increases from 0.1 to 100, the latter again agreeing with the standard result for the discontinuous DF. These results suggest that, for the standard case, the existence of a discontinuity in the DF does not play a significant role.

## 5.3 Origin of the instability islands

It is beyond the scope of this paper to characterize in detail the origin of the various manifolds. Let us nonetheless speculate on the nature of the three types of instability.

From Fig. 1, the instability island for $q < 0$ corresponds to spheres where the centrifugal support is significant and where the azimuthal pressure is low. Such systems may redistribute their orbits via an instability, in order to increase their entropy towards ‘hotter’

configurations. More generally, even for $q > 0$, mean rotation provides a reservoir of free energy for the sphere. Beyond this entropic argument, the origin of the $q < 0$ manifold is not understood in detail at this stage.

As for the $q > 0$ manifolds, the matrix method and angle-action variables allow us to emphasize the role played by orbits around specific resonances (inner Lindblad resonance (ILR), corotation (COR), and outer Lindblad resonance (OLR) for bisymmetric instabilities) in driving the instability. It provides a mathematical framework to study the tumbling instability scenario brought forward by Lynden-Bell (1979) for discs, but which also applies to spheres dominated by quasi-radial orbits. Under the assumption of the toy model presented in Pichon & Lynden-Bell (1993), which only focuses on the ILR, and assumes all orbits are identical with positive effective orbital moment of inertia<fnref n="6" /> (Earn & Lynden-Bell 1996), the instability criterion becomes the azimuthal analogue of Jeans’s (Penrose 1960). The distribution of tumbling orbits has to be narrow enough around a given maximum to allow their relative torque to overcome their spread in azimuthal tumbling rates (azimuthal pressure), which otherwise phase mixes any disturbance. The standard (non-rotating) ROI seems in most cases to qualitatively

<footnote n="6">⁶We note that the role played by the moment of inertia is implicit in equation (9) via the integral over $J$ which can be rewritten as an integral over, say, $\Omega_{ILR}$ via the corresponding moment of inertia.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Illustration of the growth rate (top panel) and oscillation frequency (bottom panel) as a function of $\alpha$ for fixed values of $q$, in the $q > 0$ region. The thin lines represent the global polynomial fits from equation (16).</figcaption></figure>

correspond to such a scenario (Polyachenko & Shukhman 2015), so long as the adiabaticity of the fast action $J_f = J_t + \frac{1}{2}L$ allows us to pin down the ILR as the unique relevant resonance (but see Polyachenko & Shukhman 2015 for a discussion of counter examples with fast growing modes).

Indeed, Fig. 14 shows that most orbits have $\zeta_{ILR} > 0$, with $\zeta_{ILR} = (\partial\Omega_{ILR}/\partial L)_{J_f}$ the cooperation parameter at the ILR, and the tumbling frequency $\Omega_{ILR} = \Omega_\phi - \frac{1}{2}\Omega_t$. Hence, the system is likely to support a tumbling instability, provided the resonance frequency distribution is sufficiently narrow (Lynden-Bell 1979). This is illustrated in Fig. 15, where the PDF of the tumbling frequency $\Omega_{ILR}$ is represented for various values of $q \ge 0$ and $\alpha$. As highlighted in these figures, the effect of increasing $q$ is both to concentrate the orbits in regions of high cooperation parameter and to lower the spread in azimuthal precession rates. Rotation also narrows the spread in azimuthal precession rates, which favours the development of instabilities. In order to take into account the direction of rotation of the orbits and make this reasoning slightly more quantitative, we classified orbits by the sign of $L_z$,<fnref n="7" /> mapping them into DF($\Omega_{ILR}$Sign($L_z$)). Without rotation ($\alpha = 0$), the DF is symmetric, since the phase space DF is independent of $L_z$. In this case, increasing $q$ concentrates the DF in the regions of low (absolute) frequency, consequently decreasing the overall spread. When $\alpha > 0$, the left bump of PDF($\Omega_{ILR}$) will be multiplied by $1 - \alpha$ and the right bump by $1 + \alpha$. This rotation-induced tilting decreases the spread in $\Omega_{ILR}$

and is likely to be responsible for the seemingly linear growth of the instability pattern speed with $\alpha$ in the Fast-ROI regime, as plotted in Fig. 10.

While increasing $q$ increases the amount of radial and cooperating orbits, decreases the spread in precession rates, and makes the system more unstable, the evolution of the growth rate with $\alpha$ is more complicated. For $\alpha > 0$, rotation provides a reservoir of rotational kinetic energy to feed possible unstable modes (Binney & Tremaine 2008). A naive expectation would be that, as one introduces positive $\alpha$ the instability gets stronger and rotates more (while preserving adiabaticity), which it does on the Slow-ROI manifold. However, for $\alpha > 0.2$, the pattern is spinning too fast to carry a large amount of orbits, which in turn lowers the destabilizing power of this process. When the pattern speed reaches the edge of the $\Omega_{ILR}$ distribution ($\alpha \simeq 0.8$ in the $q=2$ case, where $\omega_0/m \simeq 0.11$), no more orbits can contribute to the instability and it vanishes.

Orbital torquing between the COR and OLR resonances, where orbits typically have negative moment of inertia, may also make the sphere unstable, if the spread in azimuthal tumbling rates is centred on a given narrow minimum (the galactic analogue of the two-stream instability in plasma physics, Penrose 1960). In practice, the exact condition from equation (10) is more complex than Lynden-Bell's criterion since in equation (9), $\bar{M}(\omega)$ involves all couplings simultaneously.

## 6 CONCLUSION

We investigated the stability of a two-parameter sequence of modified Plummer spheres in which the velocity anisotropy (parameter $q$) and the net amount of total angular momentum (parameter $\alpha$) can be varied independently. For the first time, we implemented the matrix method to compute the growth rate, oscillation frequency, and shape of the unstable modes of a rotating stellar sphere, and compared those favourably to N-body simulations. We then charted the marginal stability boundary in the parameter space of velocity anisotropy and total angular momentum, using not only $\alpha, q$, but also $(\beta_\phi), (L_z)$, which are intended to help in applications to wider classes of model than the modified Plummer sphere.

When rotation was introduced, two apparently separate sets of modes were identified, corresponding to models with predominantly radial and tangential velocity anisotropy, respectively. For radially dominated spheres, the growing modes themselves separate into two families, which appear as two intersecting surfaces in a plot of growth rate against the two parameters, corresponding to ‘Fast’ and ‘Slow’ modes, depending on the net rotation rate of the sphere. The sequence of unstable modes that dominates much of the parameter space of tangentially biased models relies essentially on the rotational kinetic energy of the models, and may be novel. On the other hand, the $q > 0$ modes appear to extend the known ROI to rotating configurations, with a slowly rotating mode dominating at higher system rotation ($\alpha \sim 1$), perhaps surprisingly. We conjecture that the $q < 0$ modes (tangential anisotropy) correspond to the conversion of rotational support into heat, which is made possible by the low azimuthal and radial pressure. Equation (18) offers a first generalization of the classical results of Ostriker & Peebles (1973) to rotating spheres.

One possible limitation of the present mapping of instabilities is the procedure in equation (7) for the introduction of rotation, as it produces a sharp discontinuity in orbital space at $L_z = 0$. The results of Section 5.2, however, show that the discontinuity plays no significant role, at least for the standard case of $(\alpha, q) = (1, -6)$: study of a similar DF in which the discontinuity is smoothed out

<footnote n="7">⁷In such a system, all orbits precess in the same direction as they revolve.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Illustration of various modes' shapes as one explores both surfaces of instability in the $a > 0$ region. For each coloured point, the mode's density in the equatorial plane is represented. Two lines of points were chosen with constant $q$ and varying $\alpha$: $q = 0.5$, and $q = 2$. The grey arrow indicates the direction of rotation of both the global rotation and the instability patterns. On each mode, the value of $(\alpha, q)$ is written inside the circle, as well as a letter indicating whether the mode is on the Fast-ROI (F) or Slow-ROI (S) surface. In contrast to Fig. 9, the shapes of the modes of the 'Fast' and 'Slow' families are strikingly different.</figcaption></figure>

give indistinguishable results. Still, it would be of interest to study more generic classes of rotating distributions, with or without such a discontinuity, but introducing greater variety in the differential rotation and the spatial distribution of the rotation within the cluster. Such an extension of the present work is the topic of ongoing study, exploiting the flexibility of the N-body method (Breen et al. 2019b).

An additional dimension of this problem which deserves further attention concerns the detailed investigation of the role of corotation and other resonances, which, at this stage, has been only partly addressed in Section 5.3. Indeed, the behaviour of differentially

rotating polytopes (e.g. see, among others, Pickett et al. 1996; New & Shapiro 2001; Centrella et al. 2001; Shibata et al. 2002; Ou & Tohline 2006; Galeazzi et al. 2012) and the suggested connection between the emergence of an unstable mode and the existence of the associated corotation point (see e.g. Watts et al. 2005; Passamonti & Andersson 2015; Yoshida & Saijo 2017) may offer a meaningful physical intuition for the development of a comparative analysis of the stability of stellar and fluid dynamical systems.

In the present study, the focus has been on the Plummer sphere, which is appropriate to describe stellar systems with a central
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Isocontours of the growth rate in the $(\alpha, q)$-space for $q > 0$, as given by the fits from equation (16). Contours are spaced linearly between the maximum growth rate and zero. The dashed line highlights the transition between the Fast-ROI and Slow-ROI instabilities, as already visible in Fig. 11.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Illustration of the surfaces of instability in the $(\langle\beta_\theta\rangle, \langle L_z\rangle)$-space. This figure is to be compared to Fig. 7.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Superimposed isocontours of phase space DFs (equation 3) and cooperation parameter $\zeta_{ILR}$ in $(J_f, L)$ space. Filled contours: phase space DF of the $q=2$ system. The coloured regions are linearly spaced from 10 per cent (yellow) to 90 per cent (red) of the maximum value. Thin contours: phase space DF of the $q=0$ system. The contours are linearly spaced from 10 per cent (grey) to 90 per cent (black) of the maximum value. Thick contours: cooperative parameter $\zeta_{ILR}$. The contours are linearly spaced from 0 per cent (grey) to 80 per cent (black) of the maximum value. The large amount of radial cooperative orbits seems responsible for the instability of the $q=2$ systems.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** PDF of the tumbling frequency $\Omega_{ILR}$ for various values of $q \ge 0$ and $\alpha$. When $q$ increases, the PDF symmetrically concentrates around $\Omega_{ILR} = 0$. When $\alpha$ increases, some of the left hill is transferred into the right hill, concentrating the PDF around positive values of $\Omega_{ILR}$. As $\alpha$ grows, the effective frequency PDF can become narrow enough to drive a tumbling instability.</figcaption></figure>

core, but does not capture the properties of cuspy ellipticals. In that respect, the current stability charting could be extended to, e.g. rotating Hernquist models. A subsequent major step in this line of investigation could be the exploration of the properties of non-spherical equilibria (e.g. for a study of selected axisymmetric models, see Varri 2012). One difficulty in this context is the lack of angle-action variables, but Stäckel families (Stäckel 1890), with more general classes of DF than the thin shell configurations investigated in Robijn (1995), could be a possible line of attack.

In a similar vein, it would also be of interest to investigate the impact of rotation on secular evolution. The effect of such an additional
physical ingredient could be substantial, as rotation directly impacts the set of available frequencies, allowing resonant relaxation to reshuffle the system's orbital structure. The investigation presented by Hamilton et al. (2018) could be naturally extended in that respect. Certainly, N-body studies show that the effect on the time to core collapse can exceed 20 per cent (Breen et al. 2019a).

Finally, ongoing (Gaia Collaboration et al. 2018) and upcoming proper motion and spectroscopic surveys will soon provide observational constraints on the level of rotation within the several nearby globular clusters of our Milky Way. In addition, ongoing (Sydney- AAO Multi-object Integral field spectrograph, Croom et al. 2012; Mapping Nearby Galaxies at Apache Point Observatory, Bundy et al. 2015) and upcoming (Hector, Bland-Hawthorn 2015) integral field spectroscopy surveys will provide large statistical samples of the resolved kinematics of elliptical galaxies. It will, therefore, be of great interest to characterize the stability boundaries identified in the present study in terms of relevant kinematic observables, with an attention to possible degeneracies associated with projection effects. Such a development could provide a physically based framework to interpret the rich kinematics of several classes of spheroidal stellar systems which will result from upcoming empirical investigations.

## ACKNOWLEDGEMENTS

This work was motivated in part by D. Lynden-Bell's interest in the stability of rotating spheres. We dedicate this paper to his memory. JBF acknowledges support as a Hubble Fellow from Program number HST-HF2-51374, provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. This research is part of ANR grant Spin(e) (ANR-13-BS05-0005, http://cosmicorigin.org) and ERC grant 70193 (COSFORM). PGB, DCH, and ALV acknowledge support from the Leverhulme Trust (Research Project Grant, RPG- 2015-408), ALV also from a Marie Skłodowska-Curie Fellowship (MSCA-IF-EF-RI NESSY 658088). This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, cofunded by CEA and CNES. This research was supported in part by the National Science Foundation under grant no. NSF PHY-1748958. SR and ALV are grateful to the School of Mathematics and Institute for Astronomy in Edinburgh as well as the Institut d'Astrophysique de Paris for hospitality during mutual visits in Fall 2018. This work has made use of the Horizon Cluster hosted by Institut d'Astrophysique de Paris.

## REFERENCES

* Bacon R. et al., 2001, MNRAS, 326, 23

* Barnes J., Goodman J., Hut P., 1986, ApJ, 300, 112

* Bellini A., Anderson J., Bedin L. R., King I. R., van der Marel R. P., Piotto G., Cool A., 2017, ApJ, 842, 6

* Bianchini P., van der Marel R. P., del Pino A., Watkins L. L., Bellini A., Fardal M. A., Libralato M., Sills A., 2018, MNRAS, 481, 2125

* Binney J., 2014, MNRAS, 440, 787

* Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn., Princeton Univ. Press, Princeton, NJ

* Bland-Hawthorn J., 2015, in Ziegler B. L., Combes F., Dannerbauer H., Verdugo M., eds, Galaxies in 3D across the Universe, IAU Symp. Vol. 309, The Hector Survey: Integral Field Spectroscopy of 100,000 Galaxies. Kluwer, Dordrecht, p. 21

* Breen P. G. et al., 2019a, MNRAS, in prep

* Breen P. G. et al., 2019b, MNRAS, in prep

* Breen P. G., Varri A. L., Heggie D. C., 2017, MNRAS, 471, 2778

* Bundy K., Bershady M. A., Law D. R., Yan R., Drory N., MacDonald N., 2015, ApJ, 798, 7

* Cappellari M. et al., 2011, MNRAS, 413, 813

* Centrella J. M., New K. C. B., Lowe L. L., Brown J. D., 2001, ApJ, 550, L193

* Chandrasekhar S., 1969, Ellipsoidal Figures of Equilibrium. Yale Univ. Press, New Haven, CT

* Clutton-Brock M., 1973, Ap&SS, 23, 55

* Croom S. M., Lawrence J. S., Bland-Hawthorn J., Bryant J. J., Fogarty L., 2012, MNRAS, 421, 872

* Davies R., 2011, Astron. Geophys., 52, 5.18

* Davies R. L., Efstathiou G., Fall S. M., Illingworth G., Schechter P. L., 1983, ApJ, 266, 41

* Dehnen W., 2001, MNRAS, 324, 273

* Dehnen W., 2002, J. Comput. Phys., 179, 27

* Dejonghe H., 1987, MNRAS, 224, 13

* Earn D. J. D., Lynden-Bell D., 1996, MNRAS, 278, 395

* Evans N. W., An J., Bowden A., Williams A. A., 2015, MNRAS, 450, 846

* Feldmeier-Krause A., Zhu L., Neumayer N., van de Ven G., de Zeeuw P. T., Schödel R., 2017, MNRAS, 466, 4040

* Ferraro F. R. et al., 2018, ApJ, 860, 50

* Fouvry J.-B., Pichon C., Magorrian J., Chavanis P.-H., 2015, A&A, 584, A129

* Gaia Collaboration et al., 2018, A&A, 616, A1

* Galeazzi F., Yoshida S., Eriguchi Y., 2012, A&A, 541, A156

* Hamilton C., Fouvry J.-B., Binney J., Pichon C., 2018, MNRAS, 481, 2041

* Hénon M. H., 1971, Ap&SS, 14, 151

* Hernquist L., Ostriker J. P., 1992, ApJ, 386, 375

* Jeans J., 1915, MNRAS, 76, 70

* Kamann S. et al., 2018, MNRAS, 473, 5591

* Kuijken K., Dubinski J., 1994, MNRAS, 269, 13

* Lynden-Bell D., 1960, MNRAS, 120, 204

* Lynden-Bell D., 1967, MNRAS, 136, 101

* Lynden-Bell D., 1979, MNRAS, 187, 101

* Merritt D., 1999, PASP, 111, 129

* New K. C. B., Shapiro S. L., 2001, ApJ, 548, 439

* Nitadori K., Aarseth S. J., 2012, MNRAS, 424, 545

* Ostriker J. P., Peebles P. J. E., 1973, ApJ, 186, 467

* Ou S., Tohline J. E., 2006, ApJ, 651, 1068

* Passamonti A., Andersson N., 2015, MNRAS, 446, 555

* Penrose O., 1960, Phys. Fluids, 3, 258

* Pichon C., Cannon R. C., 1997, MNRAS, 291, 616

* Pichon C., Lynden-Bell D., 1993, Statistical Description of Transport in Plasma, Astro- and Nuclear Physics Les Houches Series. Nova Science Publishers, New York, p. 261

* Pickett B. K., Durisen R. H., Davis G. A., 1996, ApJ, 458, 714

* Plummer H. C., 1911, MNRAS, 71, 460

* Polyachenko V. L., Shukhman I. G., 1981, SvA, 25, 533

* Polyachenko E. V., Shukhman I. G., 2015, MNRAS, 451, 601

* Posti L., Binney J., Nipoti C., Ciotti L., 2015, MNRAS, 447, 3060

* Press W. et al., 2007, Numerical Recipes, 3rd edn., Cambridge Univ. Press, Cambridge

* Rahmati A., Jalali M. A., 2009, MNRAS, 393, 1459

* Robijn F., 1995, PhD thesis, Leiden Observatory

* Saha P., 1991, MNRAS, 248, 494

* Sanders J. L., Evans N. W., 2015, MNRAS, 454, 299

* Seth A. C., Blum R. D., Bastian N., Caldwell N., Debattista V. P., 2008, ApJ, 687, 997

* Shibata M., Karino S., Eriguchi Y., 2002, MNRAS, 334, L27

* Sollima A., Baumgardt H., Hilker M., 2019, MNRAS, 485, 1460

* Stäckel P., 1890, Math. Ann., 35, 91

* Tremaine S., 2005, ApJ, 625, 143

* Tremaine S., Weinberg M. D., 1984, MNRAS, 209, 729

* van der Marel R. P., Sigurdsson S., Hernquist L., 1997, ApJ, 487, 153

* Vandervoort P. O., 1980, ApJ, 240, 478

* Varri A. L., 2012, PhD thesis, Università degli Studi di Milano, Dipartimento di Fisica, Milano, Italy
Varri A. L., 2016, Mem. Soc. Astron. Italiana, 87, 592
Watts A. L., Andersson N., Jones D. I., 2005, ApJ, 618, L37
Weinberg M. D., 1989, MNRAS, 239, 549
Weinberg M. D., 1991, ApJ, 368, 66
Yoshida S., Saijo M., 2017, MNRAS, 466, 600

APPENDIX A: THE MATRIX METHOD

The starting point of this Appendix is equation (9), the generic expression of the response matrix of a spherically symmetric system. It involves a biorthogonal basis of potentials and densities, $\psi^{(p)}(\mathbf{x})$ and $\rho^{(p)}(\mathbf{x})$, defined as

$$ \psi^{(p)}(\mathbf{x}) = \int d\mathbf{x}' U(\mathbf{x}, \mathbf{x}') \rho^{(p)}(\mathbf{x}'), \\ \int d\mathbf{x} [\psi^{(p)}(\mathbf{x})]^* \rho^{(q)}(\mathbf{x}) = -\delta_p^q, \qquad (A1) $$

with $U(\mathbf{x}, \mathbf{x}') = -G/|\mathbf{x} - \mathbf{x}'|$ the pairwise Newtonian interaction. For a spherically symmetric system, it is natural to separate the azimuthal dependence of the basis elements with spherical harmonics and write

$$ \begin{aligned} \psi^{(p)}(\mathbf{x}) &= \psi_{lmn}(r, \vartheta, \phi) = Y_{\ell}^{m}(\vartheta, \phi) U_{n}^{\ell}(r), \\ \rho^{(p)}(\mathbf{x}) &= \rho_{lmn}(r, \vartheta, \phi) = Y_{\ell}^{m}(\vartheta, \phi) D_{n}^{\ell}(r), \end{aligned} \qquad (A2) $$

where a given basis element is characterized by three indices, $\ell \ge 0$, $|m| \le \ell$, and $n \ge 0$, and we follow the convention $\int d\vartheta d\phi \sin(\vartheta) |Y_{\ell}^{m}|^2 = 1$. We also introduced the (real) radial functions $U_{n}^{\ell}(r)$ and $D_{n}^{\ell}(r)$, and many examples can be found in the literature (Clutton-Brock 1973; Weinberg 1989; Hernquist & Ostriker 1992; Rahmati & Jalali 2009). Finally, in equation (9), for any resonance vector $\mathbf{n} \in \mathbb{Z}^3$, we introduced the Fourier transform with respect to the angle $\theta$ as

$$ f(\theta, \mathbf{J}) = \sum_{n} f_n(\mathbf{J}) e^{in\theta}, \quad f_n(\mathbf{J}) = \int \frac{d\theta}{(2\pi)^3} f(\mathbf{J}) e^{-in\theta}. \qquad (A3) $$

Following Tremaine & Weinberg (1984), the Fourier transforms of the basis elements from equation (A2) read

$$ \psi_n^{(p)}(\mathbf{J}) = \delta_{mp}^{n_3} i^{m_p-n_2} Y_{\ell_p}^{n_2}(\pi/2, 0) R_{n_2 m p}^{\ell p}(i) W_{\ell_p n p}^{\bar{n}}(\tilde{\mathbf{J}}), \qquad (A4) $$

with the reduced resonance vector $\tilde{\mathbf{n}} = (n_1, n_2)$. Here, we introduced the rotation matrix $R_{nm}^{\ell}(i)$ as

$$ R_{nm}^{\ell}(i) = \sum_{t} (-1)^t \frac{\sqrt{(\ell + n)!(\ell - n)!(\ell + m)!(\ell - m)!}}{(\ell - m - t)!(\ell + n - t)!t!(t + m - n)!} \\ \times [\cos(i/2)]^{2\ell + n - m - 2t} [\sin(i/2)]^{2t + m - n}, \qquad (A5) $$

where the sum over $t$ is restricted to the values such that the arguments of the factorials are positive, i.e. $t_{min} \le t \le t_{max}$, with $t_{min} = \text{Max}[0, n - m]$ and $t_{max} = \text{Min}[\ell - m, \ell + n]$. In equation (A4), we also introduced the Fourier-transformed 'in-plane' radial coefficients $W_{\ell n}^{\bar{n}}(\tilde{\mathbf{J}})$ defined as

$$ W_{\ell n}^{\bar{n}}(\tilde{\mathbf{J}}) = \frac{1}{\pi} \int_{0}^{\pi} d\theta_1 U_n^{\ell}(r(\theta_1)) \cos[n_1\theta_1 + n_2(\theta_2 - \psi)], \qquad (A6) $$

which are real for real radial basis functions. This integral asks for the radial dependence of the angles $\theta_1$ and $(\theta_2 - \psi)$ for a given

action $\tilde{\mathbf{J}}$. They are given by

$$ \begin{aligned} \theta_1 &= \Omega_1 \int_C dr \frac{1}{\sqrt{2(E - \psi(r)) - J_2^2/r^2}}, \\ \theta_2 - \psi &= \int_C dr \frac{\Omega_2 - J_2/r^2}{\sqrt{2(E - \psi(r)) - J_2^2/r^2}}, \end{aligned} \qquad (A7) $$

where $C$ is the integration contour going from the pericentre $r_p$ up to the current position $r = r(\theta_1)$ along the radial oscillation. Computing the radial coefficients from equation (A6) is numerically challenging. In Appendix B, we present an efficient implementation of this computation, that improves upon Fouvry et al. (2015). Finally, we note that these basis elements satisfy the symmetry relation $\psi_n^{\ell p, -m p, n p} = (-1)^{m p} (\psi_{-n}^{\ell p, m p, n p})^*.$

Let us now follow equation (9) to compute the response matrix associated with the DF from (7). Owing to its linear dependence with respect to $\alpha$, the matrices $\tilde{\mathbf{M}}_0(\omega)$ and $\tilde{\mathbf{M}}_1(\omega)$ are defined as

$$ \tilde{\mathbf{M}}_{pq}^{1}(\omega) = (2\pi)^3 \sum_{n} \int d\mathbf{J} \frac{\mathbf{n} \cdot \partial F_0 / \partial \mathbf{J}}{\omega - \mathbf{n} \cdot \Omega(\mathbf{J})} [\psi_{n}^{(p)}(\mathbf{J})]^{*} \psi_{n}^{(q)}(\mathbf{J}), \qquad (A8) $$

$$ \tilde{\mathbf{M}}_{pq}^{1}(\omega) = (2\pi)^3 \sum_{n} \int d\mathbf{J} \frac{\mathbf{n} \cdot \partial(F_0 \text{Sign}(L_z))/\partial \mathbf{J}}{\omega - \mathbf{n} \cdot \Omega(\mathbf{J})} \\ \times [\psi_{n}^{(p)}(\mathbf{J})]^{*} \psi_{n}^{(q)}(\mathbf{J}). \qquad (A9) $$

These coefficients satisfy the symmetry relations

$$ \begin{aligned} \tilde{\mathbf{M}}_{-m p, -m p}(\alpha, \omega_0 + i\eta) &= [\tilde{\mathbf{M}}_{m p, m p}(\alpha, -\omega_0 + i\eta)]^*, \\ \tilde{\mathbf{M}}_{-m p, -m p}(\alpha, \omega) &= \tilde{\mathbf{M}}_{m p, m p}(-\alpha, \omega). \end{aligned} \qquad (A10) $$

In particular, for $\alpha = 0$, we recover the invariance of the Nyquist contours with respect to complex conjugation. Indeed, an immediate consequence of equation (A10) is

$$ \tilde{\mathbf{M}}_{m p, m p}(0, \omega) = [\tilde{\mathbf{M}}_{m p, m p}(0, -\omega_0 + i\eta)]^*. \qquad (A11) $$

When a small global rotation is introduced in an originally unstable system, the Nyquist contours depart from this symmetric case, consequently shifting the values of the growth rate and oscillation frequency of instability. This behaviour is represented in Fig. A1, in the mildly unstable $q = 2, \alpha = 0$ case.

Owing to the Kronecker symbol, $\delta_{mp}^{n_3}$, in equation (A4), the sums over $n_3$ in equation (A8) can be straightforwardly carried out. In equation (A8), all the dependences with respect to $L_z$ are in the Fourier transformed basis elements, via the rotation matrices from equation (A5). We recall that they satisfy in particular the orthogonality relation

$$ \int_{0}^{\pi} d\mathbf{i} \sin(\mathbf{i}) R_{n_2 m p}^{\ell p}(\mathbf{i}) R_{n_2 m p}^{\ell q}(\mathbf{i}) = \delta_{\ell p}^{\ell q} \frac{2}{2\ell p + 1}. \qquad (A12) $$

Following the same calculation as in Hamilton et al. (2018), this allows us to immediately rewrite $\tilde{\mathbf{M}}_0(\omega)$ as

$$ \tilde{\mathbf{M}}_{pq}^{0}(\omega) = \delta_{\ell p}^{\ell q} \delta_{mp}^{m q} \sum_{\bar{\mathbf{n}}} C_{\ell p}^{n_2} P_{\ell p n p \ell p n q}^{\bar{n}}(\omega), \qquad (A13) $$

with $C_{\ell}^{n} = 2(2\pi)^3 |Y_{\ell}^{n}(\pi/2, 0)|^2/(2\ell + 1)$. In equation (A13), the coefficient $P(\omega)$ is also introduced as

$$ P_{\ell p n p \ell q n q}^{\bar{n}}(\omega) = \int d\tilde{\mathbf{J}} L \frac{\tilde{\mathbf{n}} \cdot \partial F_0 / \partial \tilde{\mathbf{J}}}{\omega - \tilde{\mathbf{n}} \cdot \tilde{\Omega}(\tilde{\mathbf{J}})} W_{\ell p n p}^{\bar{n}}(\tilde{\mathbf{J}}) W_{\ell q n q}^{\bar{n}}(\tilde{\mathbf{J}}), \qquad (A14) $$

where the prefactor $L$ originates from the volume element $d\mathbf{J} = d\tilde{\mathbf{J}} L d\mathbf{i} \sin(\mathbf{i})$. Equation (A13) corresponds to the generic
<figure><img src="image_18.png" /><figcaption>**Figure A1.** Zoom of the Nyquist contours around the instability for $q = 2$, $\alpha = 0$ (top panel) and $\alpha = 4 \times 10^{-3}$ (bottom panel). On the top panel, the Nyquist contours are symmetric with respect to the real axis, and an unstable mode with $(\eta, \omega_0) \simeq (0.0067, 0)$ is detected. On the bottom panel, rotation breaks the horizontal symmetry of the diagram. This shifts the values of growth rate and oscillation frequency to $\eta \simeq 0.009$ and $0.0045 < \omega_0 < 0.0060$.</figcaption></figure>

response matrix for a spherically symmetric and non-rotating system. One can proceed similarly to compute the matrix $\hat{M}_1(\omega)$ from equation (A9). To do so, the needed gradient can be computed as

$$ \mathbf{n} \cdot \frac{\partial(F_0(\tilde{\mathbf{J}})\mathrm{Sign}(L_z))}{\partial \mathbf{J}} = \tilde{\mathbf{n}} \cdot \frac{\partial F_0}{\partial \tilde{\mathbf{J}}} \mathrm{Sign}(L_z) + 2n_3\delta_D[L_z]F_0(\tilde{\mathbf{J}}) \quad (A15) $$

using $\partial\mathrm{Sign}(x)/\partial x = 2\delta_D(x)$, and $\mathbf{n} \cdot \partial E/\partial \mathbf{J} = \tilde{\mathbf{n}} \cdot \tilde{\Omega}(\tilde{\mathbf{J}})$. It is then straightforward to perform the integration of equation (A9). As such, we define the coefficients

$$ I_{n_2 m^p}^{\ell^p \ell^q} = (2\pi)^3 Y_{\ell^p}^{n_2} (\pi/2, 0) Y_{\ell^q}^{n_2} (\pi/2, 0) \\ \times \int_0^\pi \mathrm{d}i \sin(i) \mathrm{Sign}(\cos(i)) R_{n_2 m^p}^{\ell^p}(i) R_{n_2 m^p}^{\ell^q}(i). \quad (A16) $$

and

$$ J_{n_2 m^p}^{\ell^p \ell^q} = 2 m^p (2\pi)^3 Y_{\ell^p}^{n_2} (\pi/2, 0) Y_{\ell^q}^{n_2} (\pi/2, 0) \\ \times R_{n_2 m^p}^{\ell^p}(\pi/2) R_{n_2 m^p}^{\ell^q}(\pi/2). \quad (A17) $$

These should be seen as geometrical prefactors that can be pre-computed once and for all, and do not depend on the details of the considered system. In particular, they satisfy the symmetry

<figure><img src="image_19.png" /><figcaption>**Figure A2.** Sketch of the $\hat{M}_{pq}^0$ and $\hat{M}_{pq}^1$ matrices in the $q = -6$, $\alpha = 1$, $n_{max} = 5$, $\ell_{max} = 4$ case, with $(\eta, \omega_0) = (0.074, 0.60)$, corresponding to an $m = 2$ unstable mode. The colours scale with the norm of each cell from white, light blue up to orange. We note that $\hat{M}_{pq}^1$ is not block-diagonal anymore.</figcaption></figure>

relations

$$ \left\{ \begin{array}{l} I_{-m^p} = -I_{m^p}, \\ J_{-m^p} = -J_{m^p}, \end{array} \right. \quad \left\{ \begin{array}{l} I_{-n_2} = -I_{n_2}, \\ J_{-n_2} = J_{n_2}. \end{array} \right. \quad (A18) $$

Similarly to equation (A14), let us introduce

$$ Q_{\ell^p n^p \ell^q n^q}^{\tilde{n}}(\omega) = \int d\tilde{\mathbf{J}} \frac{F_0(\tilde{\mathbf{J}})}{\omega - \tilde{n} \cdot \tilde{\Omega}(\tilde{\mathbf{J}})} W_{\ell^p n^p}^{\tilde{n}}(\tilde{\mathbf{J}}) W_{\ell^q n^q}^{\tilde{n}}(\tilde{\mathbf{J}}), \quad (A19) $$

Gathering all these notations, the matrix $\hat{M}_1(\omega)$ reads

$$ \hat{M}_{pq}^1(\omega) = \delta_{mp}^{mq} \sum_{\tilde{n}} \left\{ P_{\ell^p n^p \ell^q n^q}(\omega) I_{n_2 m^p}^{\ell^p \ell^q} + Q_{\ell^p n^p \ell^q n^q}(\omega) J_{n_2 m^p}^{\ell^p \ell^q} \right\}. \quad (A20) $$

Figure A2 sketches $\hat{M}_{pq}^0$ and $\hat{M}_{pq}^1$ for the academic case of $n_{max} = 5$ and $\ell_{max} = 4$ (so the corresponding block is $6 \times 6$ because indexation starts at 0). Both matrices are made of blocks labelled by the pair $(\ell^p, \ell^q)$, while each element of a block is labelled by $(n^p, n^q)$. Both matrices are complex (the norm of each element is shown in Fig. A2) and symmetric, as can be seen by swapping $p$ and $q$ in equations (A13) and (A20). The off-diagonal blocks in $\hat{M}_{pq}^1$ reflects the interactions between different harmonics $\ell$. The $\delta_{\ell^p \ell^q}$ in equation (A13) makes $\hat{M}_{pq}^0$ diagonal by blocks, while the chessboard-like structure of $\hat{M}_{pq}^1$ is due to the $Y_{\ell^p}^{n_2}(\pi/2, 0) Y_{\ell^q}^{n_2}(\pi/2, 0)$ term in both $I_{n_2 m^p}^{\ell^p \ell^q}$ and $J_{n_2 m^p}^{\ell^p \ell^q}$ (see equations (A16) and (A17)), which is non-zero only when $\ell^p$ and $\ell^q$ have the same parity.

To numerically compute the integrals appearing in equation (A20), the approach of Fouvry et al. (2015) is followed, modulo two main changes: in the sampling of the orbital space $\tilde{\mathbf{J}}$, and in the computation of the elements $W_{\ell^p}(\tilde{\mathbf{J}})$, as presented in Appendix B.

In order to improve the quality of the sampling of the orbital space $\tilde{\mathbf{J}}$, we use a log-linear sampling in the $(r_p, r_a)$-space. This is done via an additional change of coordinates $(u, v) \mapsto (r_p, r_a)$, with
(u, v) sampled linearly and

$$
\begin{align}
r_p(u) &= \frac{R_b}{2(u_0 + 1)} \left[ e^{u_0 - u} \left( \operatorname{erf}\left(\frac{u_0}{\sigma}\right) + \operatorname{erf}\left(\frac{u - u_0}{\sigma}\right) \right) \right. \nonumber \\
&\qquad \left. + (u_0 - u + 1) \left( \operatorname{erfc}\left(\frac{u_0}{\sigma}\right) + \operatorname{erfc}\left(\frac{u - u_0}{\sigma}\right) \right) \right], \tag{A21} \\
r_a(u, v) &= r_p(u) + r_p(v), \nonumber
\end{align}
$$

where erf (respectively, erfc) are the (respectively, complementary) error functions, σ is an ad hoc width, and R_b is the radial scale of the basis functions. We also introduced u_0 = p_b n_max Δu - 1, with Δu the step size in the (u, v)-grid, n_max the maximum radial number for the basis functions and p_b a control parameter for the step size of the (r_p, r_a)-grid where it is linearly sampled. This transformation allows for r_p to be sampled logarithmically near 0 and linearly further out, while r_a is sampled logarithmically near r_p (close to circular orbits) and linearly further out. On top of having an explicit Jacobian, equation (A21) offers a smooth transition between the logarithmic and linear splittings, with σ = 10⁻³ the typical width of the transition region. The sampling in (r_p, r_a) becomes linear when the radial step size reaches R_b/(p_b n_max), hence the parameter p_b represents the minimal number of points in the r_p-grid per radial oscillation for the highest order radial basis function.

APPENDIX B: COMPUTING $W_{\ell n}^{\bar{n}}(\tilde{J})$

In this appendix, we detail how the radial coefficients, $W_{\ell n}^{\bar{n}}(\tilde{J})$, from equation (A6) may be efficiently computed, improving upon the method presented in Fouvry et al. (2015). We note that these coefficients only depend on the mean potential, $\psi(r)$, and our choice of basis elements, $U_n^\ell(r)$. As such, they are independent of the system's DF. They also satisfy the symmetry relation $W_{\ell n}^{\bar{n}}(\tilde{J}) = W_{\ell n}^{\bar{n}}(\tilde{J})$, which essentially halves the number of coefficients to evaluate.

In equation (A7), we have at our disposal the functions $r \mapsto \theta_1[r]$, $(\theta_2 - \psi)[r]$, so that it is more appropriate to perform the integration from equation (A6) with respect to $r$ rather than $\theta_1$. As a result, we write

$$
W_{\ell n}^{\bar{n}}(\tilde{J}) = \frac{1}{\pi} \int_{r_p}^{r_a} \frac{d\theta_1}{dr} U_n^\ell[r] \cos[n_1 \theta_1[r] + n_2(\theta_2 - \psi)[r]], \quad (B1)
$$

with the Jacobian $d\theta_1/dr = \Omega_1/\sqrt{2(E - \psi(r)) - J_2^2/r^2}$. To compute $\theta_1[r]$ and $(\theta_2 - \psi)[r]$ in equation (B1), one has to compute the integrals from equation (A7), with the boundaries $\int_C dr' = \int_{r_p}^r dr'$.

While all these integrals are well defined, we note that their integrands have a diverging behaviour at the boundaries $r = r_p$, $r_a$, originating from the vanishing of the radial velocity at the pericentre and apocentre. In order to cure these divergences and ease the integration, one can perform an additional change of variables. Rather than integrating with respect to $r$, let us therefore integrate with respect to an (explicit) angular anomaly, that we call $u$. First, we define $\Sigma = (r_p + r_a)/2$ and $\Delta = (r_a - r_p)/2$. A natural choice would be to consider the change of coordinates $r(u) = \Sigma + \Delta \sin(\pi/2 u)$, with $u \in [-1; 1]$, that is explicit, monotonic and satisfies $dr/du = 0$ for $u = \pm 1$, therefore curing the boundary divergences. In practice, it can be numerically more efficient not to use a mapping based on a trigonometric function, but rather on a polynomial function of low degree. As a result, we consider the mapping

$$
r(u) = \Sigma + \Delta f(u); \quad f(u) = u \left(3/2 - u^2/2\right), \qquad (B2)
$$

which is the appropriate third degree polynomial satisfying the four constraints $f(\pm 1) = \pm 1$ and $f'(\pm 1) = 0$ (see equation 52 in Hénon 1971). This mapping has an explicit (and simple) Jacobian and can be evaluated very efficiently. Following this change of coordinate, we rewrite equation (B1) as

$$
\begin{equation}
\begin{split}
W_{\ell n}^{\bar{n}}(\tilde{J}) = \frac{1}{\pi} \int_{-1}^{1} du \frac{dr}{du} \frac{d\theta_1}{dr} U_n^{\ell}[u] \cos[n_1 \theta_1[u]] \\
\qquad + n_2(\theta_2 - \psi)[u],
\end{split}
\tag{B3}
\end{equation}
$$

with

$$
\begin{align}
\theta_1[u] &= \int_{-1}^{u} du' \frac{dr}{du'} \frac{\Omega_1}{\sqrt{2(E - \psi[r]) - J_2^2/r^2}}, \nonumber \\
(\theta_2 - \psi)[u] &= \int_{-1}^{u} du' \frac{dr}{du'} \frac{\Omega_2 - J_2/r^2}{\sqrt{2(E - \psi[r]) - J_2^2/r^2}}. \tag{B4}
\end{align}
$$

Naively, one could interpret the expression from equation (B3) as involving nested integrals, making the numerical evaluation cumbersome. Yet, the writing of equation (B3) can be improved to make the calculation even simpler. Forgetting about indices to simplify the notations, equation (B3) asks us to compute $W(u = 1)$, for $W(u)$ of the generic form

$$
W(u) = \int_{-1}^{u} du' w[u', T_1(u'), T_2(u')], \quad (B5)
$$

where $T_1$ (respectively, $T_2$) stands for $\theta_1$ (respectively $(\theta_2 - \psi)$), and are given by expressions of the form

$$
T_1(u) = \int_{-1}^{u} du' t_1[u']; \quad T_2(u) = \int_{-1}^{u} du' t_2[u']; \qquad (B6)
$$

Because the integration boundaries of $T_1(u)$ and $T_2(u)$ match the ones for $W(u)$, one can reinterpret equation (B5) not as a set of nested integrals, but rather as only one integral. To do so, we note that $W(u)$, $T_1(u)$, and $T_2(u)$ satisfy the ordinary differential equations

$$
\begin{align}
\frac{dW}{du} &= w[u, T_1(u), T_2(u)], \nonumber \\
\frac{dT_1}{du} &= t_1[u], \nonumber \\
\frac{dT_2}{du} &= t_2[u], \tag{B7}
\end{align}
$$

with the initial conditions ($W(-1)$, $T_1(-1)$, $T_2(-1)$) = (0, 0, 0). Introducing the state vector $S = (W, T_1, T_2)$, and the velocity vector $s = (w, t_1, t_2)$, one can finally rewrite equation (B7) under the compact form

$$
\frac{dS}{du} = s[u, S(u)]; \quad S(-1) = 0. \qquad (B8)
$$

As a conclusion, under this form, computing the coefficients $W_{\ell n}^{\bar{n}}(\tilde{J})$ amounts then to computing the value of $S(1)$, by integrating forward the differential equation (B8). Following such an approach, one avoids the issue of having to compute nested integrals, and all the relevant quantities are moved forward simultaneously, which offers a significant speed-up of the numerical computation. In practice, for the applications presented in the main text, equation (B8) is integrated with a fourth-order Runge-Kutta integrator (see equation (17.1.3) in Press et al. 2007) with 100 fixed-size steps for $u \in [-1 + \varepsilon; 1 - \varepsilon]$ with $\varepsilon = 10^{-4}$ to avoid possible boundary overflows.
APPENDIX C: CONVERGENCE STUDIES

Let us briefly investigate the convergence of our numerical methods used to characterize unstable modes. Appendix C1 focuses on the

<figure><img src="image_20.png" /><figcaption>**Figure C1.** Illustration of the convergence of the matrix method for the measurement of the growth rate $\eta$ for $q = 6$ and different rotation parameters $\alpha$. Top panel: dependence of the measured growth on the maximum resonance number, $n_1^{\text{max}}$. Bottom panel: dependence of the measured growth rate with the parameter $p_b$ controlling the density of the ($r_p, r_a$)-grid, see equation (A21).</figcaption></figure>

<table><caption>**Table C1.** Growth rates and oscillation frequencies computed for various values of the parameters that control the precision of the matrix method calculations. Having chosen a particular set of parameters as a standard case, we present the computed growth rate and oscillation frequency when we vary each control parameter independently. The second column presents the interval of variation of each parameter, while the third column gives the standard value of each parameter. The last two columns show the interval of variation of the growth rate and oscillation frequency in each case. When no significant variation is measured, a single value is given. The first line gives the values of $\eta$ and $\omega_0$ for the standard case. We note that this standard case differs from the fiducial case of Section 3, which had $n_{\text{max}} = 10$ and $p_b = 2$, and is shown in the last line of the table.</caption><thead><tr><th>Parameter</th><th>Interval</th><th>Standard</th><th>η</th><th>ω<sub>0</sub></th></tr></thead><tbody><tr><td></td><td>standard</td><td></td><td>0.045</td><td>0.54</td></tr><tr><td>r<sub>min</sub></td><td>0.01–0.001</td><td>0.01</td><td>0.045</td><td>0.54</td></tr><tr><td>r<sub>max</sub></td><td>5–30</td><td>10</td><td>0.045</td><td>0.54</td></tr><tr><td>Δu</td><td>2–0.15</td><td>1</td><td>0.027–0.045</td><td>0.44–0.54</td></tr><tr><td>p<sub>b</sub></td><td>3–8</td><td>4</td><td>0.054–0.037</td><td>0.54–0.49</td></tr><tr><td>R<sub>b</sub></td><td>3–6</td><td>5</td><td>0.024–0.045</td><td>0.44–0.54</td></tr><tr><td>n<sub>1</sub><sup>max</sup></td><td>1–5</td><td>2</td><td>0.030–0.060</td><td>0.49–0.54</td></tr><tr><td>n<sub>max</sub></td><td>5–10</td><td>5</td><td>0.045–0.054</td><td>0.54–0.60</td></tr><tr><td>l<sub>max</sub></td><td>2–7</td><td>2</td><td>0.045–0.081</td><td>0.54–0.60</td></tr><tr><td></td><td>Fiducial</td><td></td><td>0.054</td><td>0.60</td></tr></tbody></table>

<figure><img src="image_21.png" /><figcaption>**Figure C2.** Illustration of the convergence of the N-body measurements for the fiducial cluster ($\alpha, q) = (1, -6$), as illustrated by the time evolution of the amplitude $|A|$ (continuous lines) and the associated exponential fit (dashed lines). Top panel: dependence with respect to the softening length $\epsilon$. This panel also includes a fully collisional run. All runs were performed with $N = 2^{17}$ particles. Bottom panel: dependence with respect to the number of particles $N$. All runs were performed with the same softening length $\epsilon = 0.01$.</figcaption></figure>

matrix method, while Appendix C2 considers the convergence of the N-body runs. In these two sections, we consider the same fiducial model as in Section 3, namely ($\alpha, q) = (1, -6$).

C1 Convergence study of the matrix method implementation

The different approximations in the computation of the response matrix introduce various control parameters that can be tuned to increase the precision (and the computational complexity) of the resulting growth rate and oscillation frequency measurements. These parameters involve on the one hand the sampling of action space in the integral of equations (A8) and (A9), which is turned into a ($r_p, r_a$)-grid by a change of variables, and on the other hand, the number of harmonics both in the Fourier decomposition in angles and in the projection onto the biorthogonal set of basis elements. The parameters controlling the ($r_p, r_a$)-grid are $r_{\text{min}}, r_{\text{max}}, \Delta u, p_b, R_b$, and were introduced in Appendix A. The parameters controlling the level of projection are $n_1^{\text{max}}$, the maximum wavenumber in the Fourier decomposition over the radial angle, $n_{\text{max}}$ the maximum order of radial basis functions, which also sets the size of a block in the response matrix, and $l_{\text{max}}$ the maximum harmonic number $l$, which also defines the number of blocks in the response matrix (see Fig. A2). We note that these two sets of parameters are not exclusive. Here, $R_b$ also controls the quality of the projection onto the radial basis set, and should be chosen so that the important information
is contained inside this radius. Similarly, $n_{max}$ also controls the grid refinement, as can be seen in equation (A21). In the following study, we made sure that varying $n_{max}$ would not change the grid sampling but only the level of projection on to the radial basis.

Table C1 presents the results of the computation of growth rates and oscillation frequencies of the fiducial model ($\alpha, q) = (1, -6)$ when we vary these control parameters. As expected (see also e.g. Hamilton et al. 2018), the results are mildly dependent on some of these parameters. Here, we show that we reached a sufficient level of convergence to confirm the existence of unstable modes, as well as the value of their growth rates and oscillation frequencies.

Figure C1 shows typical converging $\eta(\alpha)$ curves when we increase the level of precision of two of the parameters, here $p_b$ and $n_1^{max}$. In closing, our choice of control parameters ensures that for $\eta > 0.01$, the stability manifolds are properly mapped.

## C2 Convergence study of the N-body simulations

The main control parameters for the force calculation in the GYR-FALCON code are $\epsilon$ (softening parameter), $\theta$ (tolerance parameter), and the softening kernel. The softening kernel used has a density

proportional to $(1+x^2)^{-7/2}$ which, as shown by Dehnen (2001), results in much less of a force bias than Plummer softening. The parameter $\theta$ controls the accuracy of the force calculation (and is related to the opening angle in traditional tree codes). The default value of 0.6 gives relative force errors of the order of 0.001 in near-spherical systems (Dehnen 2001). Besides $N$, the softening length $\epsilon$ is the main parameter of interest since, if its value is too large, the perturbation which drives the growth of the instability might not be resolved.

In Fig. C2, we explore the effect of varying $\epsilon$ (top panel) and $N$ (bottom panel). The default value of $\epsilon = 0.05$ was found to give good agreement with results using lower values of $\epsilon$, and with a direct N-body simulation using NBODY6 (Nitadori & Aarseth 2012), which in effect has no softening, i.e. $\epsilon = 0$. The reference number of particles $N = 2^{17}$ was also found to give good agreement with larger N-body simulations: no obvious dependence of the growth rate with increasing $N$ was found.

This paper has been typeset from a TeX/LATeX file prepared by the author.