{'fraction_non_alphanumeric': 0.043590010767253655, 'fraction_numerical': 0.01722760584905005, 'mean_word_length': 4.548564270572364, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mapping the stability of stellar rotating spheres via linear response theory

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 customised 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. Generalised, approximate stability criteria are finally presented.


INTRODUCTION
Most stellar systems are rotating. While this has been clear for spiral galaxies since the beginning of extra-galactic 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 timescales, 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 w.r.t. the energy is negative, while systems with a positive gradient w.r.t. 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 (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 distribution function 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 evolution of globular clusters is now peaking, as a result of a variety of new empirical evidence. This includes recent astrometric measurements by HST (Bellini et al. 2017 (Ferraro et al. 2018) and dedicated programmes with MUSE (Kamann et al. 2018). Similarly, line-of-sight measurement 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 /σ vs. 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 et al. 1996;New & Shapiro 2001;Centrella et al. 2001;Shibata et al. 2002;Ou & Tohline 2006 and, more recently, Galeazzi et al. 2012; suggestions for a possible dynamical interpretation in which corotation points play a primary role have been offered, among others, by Watts et al. 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 et al. 2017Breen et al. , 2019a. 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., see Barnes et al. 1986;Kuijken & Dubinski 1994;van der Marel et al. 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.
The present paper is organised as follows. Section 2 presents the equilibrium model which we consider as initial state, in particular its distribution function (DF) with velocity anisotropies, and how solid rotation is introduced (through the "Lynden-Bell trick"). With the aim of characterising 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 Nbody simulations. Section 4 extends this analysis to the whole twodimensional parameter space of velocity anisotropy and rotation, and investigates their respective roles in sourcing dynamical insta-bilities. We finally discuss our results in Section 5, and conclude in Section 6.

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 distribution function that encodes both the velocity anisotropies and the system's total rotation.

Angle-action coordinates
The mean stellar system is always assumed to be spherically symmetric, and characterised by a mean radial potential ψ(r), with (r, ϑ, φ) the spherical coordinates. Throughout this study, we take the mean potential to be a Plummer potential (Plummer 1911), and work in the dimensionless units G = Rs = Mtot = 1, with Mtot the system's total mass. A spherically symmetric potential guarantees the existence of angle-action coordinates, (θ, J), (Binney & Tremaine 2008), with 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) = Lz/L, with 0 i π.
According to the Jeans theorem (Jeans 1915), the mean-field distribution of the system is fully characterised by its DF, F (J), normalised as dθdJ F = Mtot. Because the system is spherically symmetric, the orbital frequencies, Ω(J) = ∂H/∂J, satisfy the degeneracy condition Ω3(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 J = (Jr, L), and the associated reduced frequencies Ω = ∂H( J)/∂ J. Note that spherical systems with a spherically symmetric velocity distribution satisfy F = F ( J).

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) with q controlling the flavour and degree of velocity anisotropy of the system, and Γ the gamma function. Here, F is expressed in terms 1 of the hypergeometric function 2F1 , if x 1. (4) As illustrated in Figure 1, the parameter q controls the flavour of the velocity anisotropy in the system, with q < 0 (resp. q > 0) associated with tangential (resp. radial) velocity anisotropies. One of 10 -2 10 -1 10 0 10 1 r 0 α=0.5 α=1 Figure 1. Illustration of the radial dependence of v φ (r) (top panel) and β φ (r) (bottom panel) for various velocity anisotropies (via q, indicated by the shade of colour) and total angular momentum (via α, indicated by the line style). Here, α controls v φ (r) , while q controls β φ , noticeably in the outer region of the cluster. We note that for (α, q) = (1, −6) a significant fraction of the radial support is centrifugal, whereas σ φ remains comparatively small. The sphere has therefore a large reservoir of kinetic energy, which can be converted into thermal energy through an instability. the advantages of this DF is the simple forms taken by its velocity dispersions in all directions. Indeed, one has σ 2 r (r) = (1 + r 2 ) −1/2 6 − q , leading to the anisotropy parameter: We note that after introducing rotation as below, this relation still holds, as long as one uses equation (6) to define β θ , ignoring the dispersion in φ.

Introducing rotation
Following Lynden-Bell (1962), a simple way of introducing rotation in a spherically symmetric system is to modify its DF as 2 F (α, q, J) = F0(q, J) + αF0(q, J) Sign(Lz) , 2 We note that any odd fonction of Lz could be used to generalise our choice of Sign(Lz) to introduce rotation. Figure 2. Illustration of Lynden-Bell's trick with two orbits, red and blue, having the same total angular momentum modulus (the red and blue arrows are perpendicular to the orbital planes), but opposite Lz (black arrows). The α parameter of equation (7) defines the relative number of such orbits (which only differ by the orientation of their orbital plane). When α = 1, no blue orbits are left.
with the sign function Sign(x) equal to −1 for x < 0, then 0 for x = 0 and 1 for x > 0. In equation (7), |α| 1 is a dimensionless parameter that controls the amount of rotation in the system. We note that F (α = 0, q, J) = F0(q, J), corresponding to the case without rotation, and that the normalization of F is guaranteed by that of F0, since L −L dLz Sign(Lz) = 0. The DF of equation (7) does not change the distribution of the orbits' shapes, but rather the distribution of their orientations, as illustrated in Figure 2. For each orbit (i.e. each J), the probability to have Lz > 0 is multiplied by 1 + α, compared to the non-rotating DF. In the limiting case α = 1, all particles have Lz > 0.
Since a more extensive characterisation of the properties of these configurations will be provided in Breen et al. (2019a), here we just wish to highlight a few selected features 3 . In Figure 1, we illustrate how α creates some mean azimuthal motion in an otherwise non-rotating system. We also note that α induces an anisotropy in the azimuthal motions, by breaking the symmetry in the two tangential directions, as given by the azimuthal anisotropy parameter The bottom panel of Figure 1 also shows how α 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 = 0, a radial variation of β φ is also introduced in the outer region of the cluster. Moreover, increasing q decreases β φ . We note that by construction v φ (r) ∝ α. Finally, we note that even for q = −6, when α = 1, σr > σ φ . In Figure 3, we illustrate the PDF of (vr, v φ ) for various values of α and q. We note that the ansatz given by equation (7) introduces a sharp discontinuity in the distribution of stars in the v φ -direction as soon as α = 0. It is possible that this discontinuity plays a role in driving the instabilities, as discussed further in Section 5.2. . Number density of stars in the (vr, v φ )-space at the constant radius r = 1. Whenever α = 0, the PDF is discontinuous on the line v φ = 0.
As q increases, the spread in v φ decreases, and the system becomes more radially biased. As α increases, stars are swapped towards positive v φ .

Matrix method and rotating spheres
The stability properties of anisotropic, spherical equilibria have already been the object of many studies (Polyachenko & Shukhman 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 characterised by its response matrix M(ω) (Binney & Tremaine 2008), which in the case of 3D spherical systems is given by A system is said to be linearly unstable if there exists a frequency ω = ω0 + iη (with η > 0), such that M(ω) admits an eigenvalue equal to 1. The instability criterion for the sphere therefore reads In that case, the system supports an unstable mode of oscillation frequency ω0 and growth rate η. To capture the system's selfgravitating response, equation (9) involves a biorthogonal basis of potentials and densities, ψ (p) (x) and ρ (p) (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 where M0(ω) corresponds to the system's response in the absence of rotation, while α M1(ω) captures the effects associated with rotation. Equation (11) is the master equation characterising its linear stability 4 . Because of its reduced number of symmetries, we also note that the matrix M1(ω) is less sparse than the zero-rotation matrix M0(ω), and hence the numerical task is more cumbersome. In Appendix B, we detail our numerical implementation of the computation of M(α, ω). Before applying the previous machinery to the study of the linear stability of a wide range of equilibria (i.e. different (α, q)), we will start by considering one particular case. Following equation (7), we focus on the fiducial model (α, q) = (1, −6), i.e. a model with the maximum rotation and with tangential velocity anisotropies, as already illustrated in Figure 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 bi-symmetric modes dominate.

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 η > 0, one constructs the curve ω0 → det[I − M(ω0 + iη)]. If ever this contour encloses the origin, then the cluster supports an unstable mode with a growth rate larger than η. Figure 4 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 w.r.t. these various numerical parameters. In order to find the instability with the largest growth rate, we draw Nyquist contours for progressively decreasing values of η, 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 (ω0, η) (0.547, 0.047). We note that as one keeps decreasing the value of η, additional sub-dominant 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 with Xp the eigenvector. In Figure 5, we represent in the equatorial plane the density of the mode of the fiducial model. In Figure 6, 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.

N -body methodology
In order to validate the linear stability analysis further, Nbody simulations were carried out for the same fiducial model, values 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 M 0 (ω), the shift in eigenvalues and eigenvectors associated with the introduction of α 1 can be computed perturbatively. The top panel represents the Nyquist diagram without renormalization. Following Pichon & Cannon (1997), in the bottom panel, a logarithmic scaling is applied, re iθ → 1 3 log(1+10 3 r) e iθ , 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 ω 0 along each curve are indicated with the gray dashed lines. This diagram shows that the cluster supports an unstable mode with a growth rate, 0.045 η 0.049, and oscillation frequency, 0.547 ω 0 0.550.
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 characterise the frequencies of the unstable modes, we Fourier-analysed the space density as follows. For a given pattern number m = 2, we define with ρ(R, z, φ, 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 numericallymeasured function t → Am(t), as one has In practice, measuring the growth rate of an instability in Nbody data has two major difficulties. The first is noise due to the dis-  Isocontours are logarithmically sampled, with a factor ∼ 1.6 between two consecutive lines, so that the outermost line has a density 1.5×10 4 times smaller than the peak value. The similarity between the two results is striking.
crete 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 distribution function. 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.
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 |Am(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 normalised amplitude) falls below a predetermined value of tolerance. The choice of tolerance 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 that the interval is not reduced merely by the presence of noise. Based on a grid of 300 different realisations of the reference model (N = 2 17 , q = −6, α = 1), a value of 0.009 was found to produce robust measurements of the growth rate η. Similarly, ω0 is determined by a linear fit to Arg(Am(t)) (mod 2π to account for complete revolutions) over the same time interval. Figure 7 illustrates such measurements of η and ω0 for the fiducial model.
Having characterised 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 onto the same ρ (p) (r, ϑ, φ) basis as for the matrix method, to reconstruct the mode using equation (12). This is illustrated in the bottom panel of Figure 5. In Appendix C2, we briefly check the convergence of the N -body runs w.r.t. their various control parameters.
The comparison between linear theory and N -body simulations in the fiducial model shows a ∼ 15% discrepancy in the measured value of η, and a ∼ 3% discrepancy in ω0. Moreover, as highlighted in Figure 5, 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. Arg(Am(t)) (right panel) for the same cluster model as in Figure 4. Following equation (14), linear fits are also represented. We recovered (ω 0 , η) = (0.57, 0.022) that should be compared with the measurements of Figure 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 (α, 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 η, the oscillation frequency ω0, and the mode's shape.

Instability in the (α, q)-space
We reproduced the methodology of Section 3.2 to chart the (α, q)space using the linear matrix method, as well as N -body simulations. The resulting frequencies are illustrated in Figure 8.
In that figure, one can straightforwardly identify two types of instabilities. The first one, corresponding to the pyramid on the top of the η(α, q) panel (q < 0, α 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 be also be noted in the bottom right part of the η(α, 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 panels of Figure 8 show the N -body counterparts of the measurements using the matrix method. Again, the quantitative matchings between the left and right panels validate our generalisation of the matrix method. On all these panels, for both η and ω0, the lower blue plane corresponds to a growth rate threshold, η = 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) 5 .
Having validated the methodology against N -body simulations, we may now carry out a more thorough exploration of the (α, q)-space using only the matrix method, successively for tangentially-biased (q < 0) and radially-biased systems (q > 0).
(a) Growth rate, matrix method.
(c) Oscillation frequency, matrix method. (d) Oscillation frequency, N -body measurements. Figure 8. Illustration of the dependence of the growth rates η (top panels) and oscillation frequencies ω 0 (bottom panels), as a function of the cluster's parameters (α, q), measured using the response matrix method (left panels) and N -body methods (right panels). We searched for unstable modes on a (α, q)-grid composed of the locations α = 0, 0.2, 0.4, 0.6, 0.8, 1 and q = −16, −12, −6, −2, 0, 1, 2. The blue plane represents our growth rate threshold, η = 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, α > 1 2 ) and an extension of the known ROI instability island on the right (near q ∼ 2). The comparisons of the matrix method's and N -body measurements highlights again the reliability of our implementation of the linear theory.
4.2 Instability of tangentially-biased systems: q < 0 In Figure 9, we illustrate cuts of constant q in the η(α, q) and ω(α, q) surfaces in the q < 0 region. In the relevant region, the instability surface (α, q) → ω0, η can be fitted globally by the simple polynomial expressions whose cuts are also represented in Figure 9.
In Figure 10, we represent the η(α, q) surface for tangentiallybiased systems. We emphasise that in this regime, both η and ω0 increase with α, and increase as q decreases. We also note that, so far, these instabilities appear to emerge only in the presence of nonvanishing angular momentum (i.e., they are not present in configuration with α → 0). We also emphasise that their specific behaviour depend on the adopted choice of the angular momentum distribution (i.e., as resulting from the application of the "Lynden-Bell trick"), as discussed in Section 5.2. In Figure 10, 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 Figure 11, we illustrate cuts of constant q in the η(α, q) and ω(α, q) surfaces in the q > 0 region. Similarly, in Figure 12, we explore the η(α, 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 (α 1 2 ). It is represented with the solid lines in Figure 11, and the dark grey surface in Figure 12. It can be identified as the continuation in the rotating regime of the standard radial orbit instability (ROI, Polyachenko & Shukhman 1981). This type will be referred to hereafter as "Fast-ROI", since the 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 Figure 12.
The second type of instability dominates for fast-rotating systems (α 1 2 ). It is represented with the dashed lines in Figure 11, and the light grey surface in Figure 12. We note that the oscillation frequencies of these modes are small, and we therefore refer to this instability as "Slow-ROI". As illustrated in Figure 12, 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 (α, q) → ω0, η can be fitted globally by simple polynomial functions as for the Fast-and the Slow-ROI respectively. In Figure 13, we finally illustrate isocontours of the growth rates of the dominant instabilities in the (α, q > 0)-space. This fig-ure clearly highlights the transition between the two regions of instability, Fast-ROI vs. Slow-ROI.

DISCUSSION
The previous sections identified various islands of instability in the two parameters (α, 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 parametrisation 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.

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 β θ and Lz respectively, using equations (6) and (7). They read 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 α are direct proxy for β θ and Lz . Using this reparametrisation, we present in Figure 14 a remapping of the stability manifolds in the ( β θ , Lz )-space.
As expected the qualitative differences between Figures 8  and 14 are minor. The following simple fits represent the three manifolds respectively for negative q, Fast-ROI, Slow-ROI, using the shortening notations β = β θ and L = Lz . From equation (18), the Ostriker-Peebles criterion (Ostriker & Peebles 1973) can be generalised by requesting that η is above a given threshold ηcrit. For instance, for the Slow-ROI instability, it reads Lz /7 + β θ /10 ηcrit + 1 21 .

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 the Lynden-Bell trick from equation (7) which may partially impact artificially the distribution of orbits near Lz = 0, as can be seen in Figure 3. This may prevent other instabilities from playing a role within the top left part of the chart from Figure 8, though it should be noted that clusters with q ∼ 0 and α ∼ 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, Q was set to zero in equation (A19), effectively eliminating the contribution of the discontinuity of Sign(Lz) in equation (7). All three island are still present, though the growth rates are lower. As a further check, a distribution whose odd part in Lz scales like Lz/L was also investigated. Once again, qualitatively, the three instability islands remain.

Origin of the instability islands
It is beyond the scope of this paper to characterise in detail the origin of the various manifolds. Let us nonetheless speculate on the nature of the three types of instability.
From Figure 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 emphasise the role played by orbits around specific resonances (ILR, COR and OLR for bi-symmetric 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 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 correspond to such a scenario (Polyachenko & Shukhman 2015), so long as the adiabaticity of the fast action J f = Jr + 1 2 L allows us to pin down the ILR as the unique relevant resonance (but Indeed, Figure 15 shows that most orbits have ζILR > 0, with ζILR = (∂ΩILR/∂L)J f the cooperation parameter at the ILR, and the tumbling frequency ΩILR = Ω φ − 1 2 Ωr. 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 Figure 16, where the PDF of the tumbling frequency ΩILR is represented for various values of q 0 and α. 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 Lz 7 , mapping them into DF(ΩILRSign(Lz)). Without rotation (α = 0), the DF is symmetric, since the phase space DF is independent of Lz. In this case, increasing q concentrates the DF in the regions of low (absolute) frequency, consequently decreasing the overall spread. When α > 0, the left bump of PDF(ΩILR) will be multiplied by 1−α and the right bump by 1+α. This rotation-induced tilting decreases the spread in ΩILR and is likely to be responsible for the seemingly linear growth of the instability pattern speed with α in the Fast-ROI regime, as plotted in Fig. 11.
While increasing q increases the amount of radial and coop-7 In such a system, all orbits precess in the same direction as they revolve.
erating orbits, decreases the spread in precession rates, and makes the system more unstable, the evolution of the growth rate with α is more complicated. For α > 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 α the instability gets stronger and rotates more (while preserving adiabaticity), which it does on the Slow-ROI manifold. However, for α > 0.2, the pattern is spinning too fast to carry a large amount of orbits, which in turn lowers the destabilising power of this process. When the pattern speed reaches the edge of the ΩILR distribution (α 0.8 in the q = 2 case, where ω0/m 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), M(ω) involves all couplings simultaneously.

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 α) 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 α, q, but also β θ , Lz , 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 (α ∼ 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 generalisation 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 Lz = 0. It would be of interest to study a more generic class of rotating distributions, without such a discontinuity (see, for example, Section 5.2). The lack of control of differential rotation at fixed mean rotation is also a limitation. An extension of the present work to a more general class of models parametrising the shear and smoothing the discon- Figure 12. Illustration of various modes' shapes as one explores both surfaces of instability in the q > 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 α: 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 (α, 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 Figure 10, the shapes of the modes of the "Fast" and "Slow" families are strikingly different. tinuity is the topic of ongoing work, 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 polytropes (e.g., see, among others, Pickett et al. 1996 In the present study the focus has been on the Plummer sphere, which is appropriate to describe stellar systems with a central core, but it does not capture properties of cuspy ellipticals. In that re- Figure 13. Isocontours of the growth rate in the (α, 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 Figure 12. spect, 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 nonspherical 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% (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 (SAMI, Croom et al. 2012, MANGA, 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 characterise 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.
JBF acknowledges support as a Hubble Fellow from Program number HST-HF2-51374 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the As-   (3)) and cooperation parameter ζ ILR in (Jr, L) space. Filled contours: phase space DF of the q = 2 system. The coloured regions are linearly spaced from 10% (yellow) to 90% (red) of the maximum value. Thin contours: phase space DF of the q = 0 system. The contours are linearly spaced from 10% (grey) to 90% (black) of the maximum value. Thick contours: Cooperative parameter ζ ILR . The contours are linearly spaced from 0% (grey) to 80% (black) of the maximum value. The large amount of radial cooperative orbits seems responsible for the instability of the q = 2 systems.  Figure 16. PDF of the tumbling frequency Ω ILR for various values of q 0 and α. When q increases, the PDF symmetrically concentrates around Ω ILR = 0. When α increases, some of the left hill is transferred into the right hill, concentrating the PDF around positive values of Ω ILR . As α grows, the effective frequency PDF can become narrow enough to drive a tumbling instability.
where the sum over t is restricted to the values such that the arguments of the factorials are positive, i.e. tmin t tmax, with tmin = Max[0, n − m] and tmax = Min[ − m, + n]. In equation (A4), we also introduced the Fourier-transformed "in-plane" radial coefficients W n n ( J) defined as which are real for real radial basis functions. This integral asks for the radial dependence of the angles θ1 and (θ2 − ψ) for a given action J. They are given by where C is the integration contour going from the pericentre rp up to the current position r = r(θ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 ψ Let us now follow equation (9) to compute the response matrix associated with the DF from (7). Owing to its linear dependence w.r.t. α, the matrices M0(ω) and M1(ω) are defined as These coefficients satisfy the symmetry relations In particular, for α = 0, we recover the invariance of the Nyquist contours w.r.t. complex conjugation. Indeed, an immediate consequence of equation (A9) is 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 Figure A1, in the mildly unstable q = 2, α = 0 case. Owing to the Kronecker symbol, δ n 3 m p , in equation (A4), the sums over n3 in equation (A8) can be straightforwardly carried out. In equation (A8), all the dependences w.r.t. Lz are in the Fourier transformed basis elements, via the rotation matrices from equation (A5). We recall that they satisfy in particular the orthogonality relation Following the same calculation as in Hamilton et al. (2018), this allows us to immediately rewrite M0(ω) as with C n = 2 (2π) 3 |Y n (π/2, 0)| 2 /(2 + 1). In equation (A12), the coefficient P (ω) is also introduced as P n p n p q n q (ω) = d JL where the prefactor L originates from the volume element dJ = d JLdi sin(i). Equation (A12) corresponds to the generic response matrix for a spherically symmetric and non-rotating system. One can proceed similarly to compute the matrix M1(ω) from equation (A8). To do so, the needed gradient can be computed as Figure A2. Sketch of the M 0 pq and M 1 pq matrices in the q = −6, α = 1, nmax = 5, max = 4 case, with (η, ω 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 M 1 pq is not block-diagonal anymore.
(A14) using ∂Sign(x)/∂x = 2δD(x), and n · ∂E/∂J = n · Ω( J). It is then straighforward to perform the integration of equation (A14). As such, we define the coefficients Y n 2 p (π/2, 0) Y n 2 q (π/2, 0) term in both I p q n 2 m p and J p q n 2 m p (equations (A15) and (A16)), which is non-zero only when p and q have the same parity.
To numerically compute the integrals appearing in equation (A19), the approach of Fouvry et al. (2015) is followed, modulo two main changes: in the sampling of the orbital space J, and in the computation of the elements W n n ( J), as presented in Appendix B.
In order to improve the quality of the sampling of the orbital space J, we use a log-linear sampling in the (rp, ra)−space. This is done via an additional change of coordinates (u, v) → (rp, ra), with (u, v) sampled linearly and where erf (resp. erfc) are the (resp. complementary) error functions, σ is an ad hoc width, R b is the radial scale of the basis functions. We also introduced u0 = p b nmax ∆u − 1, with ∆u the step size in the (u, v)-grid, nmax the maximum radial number for the basis functions and p b a control parameter for the step size of the (rp, ra)-grid where it is linearly sampled. This transformation allows for rp to be sampled logarithmically near 0 and linearly further out, while ra is sampled logarithmically near rp (close to circular orbits) and linearly further out. On top of having an explicit Jacobian, equation (A20) offers a smooth transition between the logarithmic and linear samplings, with σ = 10 −3 the typical width of the transition region. The sampling in (rp, ra) becomes linear when the radial step size reaches R b /(p b nmax), hence the parameter p b represents the minimal number of points in the rp-grid per radial oscillation for the highest order radial basis function.
APPENDIX B: COMPUTING W n n ( J) In this Appendix, we detail how the radial coefficients, W n n ( 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, ψ(r), and our choice of basis elements, U n (r). As such, they are independent of the system's DF. They also satisfy the symmetry relation W − n n ( J) = W n n (J), which essentially halves the number of coefficients to evaluate.
In equation (A7), we have at our disposal the functions r → θ1[r], (θ2 − ψ) [r], so that it is more appropriate to perform the integration from equation (A6) w.r.t. r rather than θ1. As a result, we write with the Jacobian dθ1/dr = Ω1/ 2(E −ψ(r))−J 2 2 /r 2 . To compute θ1[r] and (θ2 −ψ)[r] in equation (B1), one has to compute the integrals from equation (A7), with the boundaries C dr = r rp dr . While all these integrals are well-defined, we note that their integrands have a diverging behaviour at the boundaries r = rp, ra, 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 w.r.t. r, let us therefore integrate w.r.t. an (explicit) angular anomaly, that we call u. First, we define Σ = (rp +ra)/2 and ∆ = (ra −rp)/2. A natural choice would be to consider the change of coordinates r(u) = Σ + ∆ sin( π 2 u), with u ∈ [−1; 1], that is explicit, monotonic and satisfies dr/du = 0 for u = ±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 which is the appropriate third degree polynomial satisfying the four constraints f (±1) = ±1 and f (±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 with 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 where T1 (resp. T2) stands for θ1 (resp. (θ2 − ψ)), and are given by expressions of the form Because the integration boundaries of T1(u) and T2(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), T1(u), and T2(u) satisfy the ordinary differential equations with the initial conditions (W (−1), T1(−1), T2(−1)) = (0, 0, 0). Introducing the state vector S = (W, T1, T2), and the velocity vector s = (w, t1, t2), one can finally rewrite equation (B7) under the compact form dS du = s u, S(u) ; S(−1) = 0.
As a conclusion, under this form, computing the coefficients W n n ( 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 fixedsize steps for u ∈ [−1 + ε; 1 − ε] with ε = 10 −4 to avoid possible boundary overflows.

APPENDIX C: CONVERGENCE STUDIES
Let us briefly investigate the convergence of our numerical methods used to characterise unstable modes. Section C1 focuses on the matrix method, while Section C2 considers the convergence of the N -body runs. In these two sections, we consider the same fiducial model as in Section 3, namely (α, 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 equation (A8), which is turned into a (rp, ra)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 (rp, ra)-grid are rmin, rmax, ∆u, p b , R b and were introduced in Appendix A. The parameters controlling the level of projection are n max 1 , the maximum wave number in the Fourier decomposition over the radial angle, nmax the maximum order of radial basis functions, which also sets the size of a block in the response matrix, and max the maximum harmonic number , which also defines the number of blocks in the response matrix (see Figure 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, nmax also controls the grid refinement, as can be seen in equation (A20). In the following study, we made sure that varying nmax would not change the grid sampling but only the level of projection onto the radial basis. Table C1 presents the results of the computation of growth rates and oscillation frequencies of the fiducial model (α, 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 η(α) curves when we increase the level of precision of two of the parameters, here p b and n max 1 . In closing, our choice of control parameters ensures that for η > 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 gyrfalcON code are ε (softening parameter), θ (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  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 η and ω 0 for the standard case. We note that this standard case differs from the fiducial case of Section 3, which had nmax = 10 and p b = 2, and is shown in the last line of the    Figure C2. Illustration of the convergence of the N -body measurements for the fiducial cluster (α, 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 w.r.t. the softening length ε. This panel also includes a fully collisional run. All runs were performed with N = 2 17 particles. Bottom panel: Dependence w.r.t. the number of particles N . All runs were performed with the same softening length ε = 0.01.
(2001), results in much less of a force bias than Plummer softening. The parameter θ 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 ε is the main parameter of interest since, if its value is too large, the perturbation which drives the growth of instability might not be resolved.
In Figure C2, we explore the effect of varying ε (top panel) and N (bottom panel). The default value of ε = 0.05 was found to give good agreement with results using lower values of ε, and with a direct N -body simulation using NBODY6 (Nitadori & Aarseth 2012), which in effect has no softening, i.e. ε = 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.