Mapping the stability of stellar rotating spheres via linear response theory

TL;DR

This paper investigates the stability of rotating spherical stellar systems using linear response theory, revealing how rotation and anisotropy influence unstable modes and providing stability criteria.

Contribution

It adapts the response matrix method to spherical rotating systems and maps stability boundaries considering anisotropy and rotation.

Findings

Identification of unstable modes related to anisotropy and rotation.

Mapping of marginal stability boundaries in parameter space.

Comparison of theoretical predictions with N-body simulations.

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

Figures

Tables

Table 1. Table 1: 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 ω 0 subscript 𝜔 0 \omega_{0} for the standard case. We note that this standard case differs from the fiducial case of Section  3 , which had n max = 10 subscript 𝑛 max 10 {n_{\mathrm{max}}=10} and p b = 2 subscript 𝑝 b 2 {p_{\mathrm{b}}=2} , and is shown in the last line of the table.

Parameter	Interval	Standard	$\eta$	${\omega }_0$
	Standard		0.045	0.54
$r_{\min }$	0.01 - 0.001	0.01	0.045	0.54
$r_{\max }$	5 - 30	10	0.045	0.54
$\Delta u$	2 - 0.15	1	0.027 - 0.045	0.44 - 0.54
$p_b$	3 - 8	4	0.054 - 0.037	0.54 - 0.49
$R_b$	3 - 6	5	0.024 - 0.045	0.44 - 0.54
$n_1^{\max }$	1 - 5	2	0.030 - 0.060	0.49 - 0.54
$n_{\max }$	5 - 10	5	0.045 - 0.054	0.54 - 0.60
${\ell }_{\max }$	2 - 7	2	0.045 - 0.081	0.54 - 0.60
	Fiducial		0.054	0.60

Equations

$$\psi (r)=-\frac{G{M}_{\mathrm{tot}}}{R_{\mathrm{s}}}\frac{1}{\sqrt{1+(r/R_{\mathrm{s}}{)}^2}},$$

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

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

$$\mathcal{F}(a,b,c,d;x)=\{\begin{array}{rlr}& x^a\frac{{}_2{F}_1(a+b,1+a-c;a+d;x)}{\Gamma (c-a)\Gamma (a+d)},& \text{if }x\le 1;\\ & \frac{1}{x^b}\frac{{}_2{F}_1(a+b,1+b-d;b+c;\frac{1}{x})}{\Gamma (d-b)\Gamma (b+c)},& \text{if }x\ge 1.\end{array}$$

$${\sigma }_r^2(r)$$

$${\sigma }_{\theta }^2(r)$$

$${\beta }_{\theta }(r)=1-\frac{{\sigma }_{\theta }^2(r)}{{\sigma }_r^2(r)}=\frac{q}{2}\frac{r^2}{1+r^2}.$$

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

$${\beta }_{\varphi }(r)=1-\frac{{\sigma }_{\varphi ,\alpha }^2(r)}{{\sigma }_r^2(r)}.$$

$$\!\widehat{M}_{pq}(\omega)\!=\!(2\pi)^{3}\!\!\sum_{\mathbf{n}\in\mathbb{Z}^{3}\,}\!\!\int\!\!\mathrm{d}\mathbf{J}\,\frac{\mathbf{n}\cdot\partial F/\partial\mathbf{J}}{\omega-\mathbf{n}\cdot\mathbf{\Omega}(\mathbf{J})}\big{[}\psi^{(p)}_{\mathbf{n}}(\mathbf{J})\big{]}^{*}\!\psi_{\mathbf{n}}^{(q)}(\mathbf{J}).$$

$$\det [\mathbf{I}-\hat{\mathbf{M}}({\omega }_0+\mathrm{i}\eta )]=0.$$

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

$$\rho (r,\vartheta ,\varphi )=\mathrm{Re}\left(\sum _p{X}_p{\rho }^{(p)}(r,\vartheta ,\varphi )\right),$$

$$C_m(t)={\int }_0^{+\infty }\mathrm{d}rr{\int }_{-\infty }^{+\infty }\mathrm{d}z{\int }_0^{2\pi }\frac{\mathrm{d}\varphi }{2\pi }\rho (r,z,\varphi ,t){\mathrm{e}}^{-\mathrm{i}m\varphi },$$

$$\frac{\mathrm{d}\log (|A_m(t)|)}{\mathrm{d}t}=\eta ;\frac{\mathrm{d}\mathrm{Arg}(A_m(t))}{\mathrm{d}t}={\omega }_0.$$

$${\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}.$$

$${\omega }_0=\frac{4\alpha }{11}-\frac{q}{34}+\frac{1}{38},$$

$${\omega }_0=\frac{2\alpha }{25}-\frac{q}{123}+\frac{1}{39},$$

$$⟨{\beta }_{\theta }⟩$$

$$⟨L_z⟩$$

$${\omega }_0=\frac{10L}{7}-\frac{\beta }{43}-\frac{1}{7},$$

$${\omega }_0=L,$$

$${\omega }_0=\frac{L}{3},$$

$${\psi }^{(p)}(\mathbf{x})=\int \mathrm{d}{\mathbf{x}}^{\prime }U(\mathbf{x},{\mathbf{x}}^{\prime }){\rho }^{(p)}({\mathbf{x}}^{\prime }),$$

$$\displaystyle\,\!\!\int\!\!\mathrm{d}\mathbf{x}\,\big{[}\psi^{(p)}(\mathbf{x})\big{]}^{*}\,\rho^{(q)}(\mathbf{x})=-\delta_{p}^{q},$$

$${\psi }^{(p)}(\mathbf{x})={\psi }_{\ell mn}(r,\vartheta ,\varphi )=Y_{\ell }^m(\vartheta ,\varphi )U_n^{\ell }(r),$$

$${\rho }^{(p)}(\mathbf{x})={\rho }_{\ell mn}(r,\vartheta ,\varphi )=Y_{\ell }^m(\vartheta ,\varphi )D_n^{\ell }(r),$$

$$f(\mathbf{\theta },\mathbf{J})=\sum _{\mathbf{n}}{f}_{\mathbf{n}}(\mathbf{J}){\mathrm{e}}^{\mathrm{i}\mathbf{n}\cdot \mathbf{\theta }},f_{\mathbf{n}}(\mathbf{J})=\int \frac{\mathrm{d}\mathbf{\theta }}{(2\pi {)}^3}f(\mathbf{J}){\mathrm{e}}^{-\mathrm{i}\mathbf{n}\cdot \mathbf{\theta }}.$$

$${\psi }_{\mathbf{n}}^{(p)}(\mathbf{J})={\delta }_{m^p}^{n_3}{\mathrm{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}^{\tilde{\mathbf{n}}}(\tilde{\mathbf{J}}),$$

$$R_{nm}^{\ell }(i)=$$

$$\displaystyle\;\;\times\big{[}\cos(i/2)\big{]}^{2\ell+n-m-2t}\,\big{[}\sin(i/2)\big{]}^{2t+m-n},$$

$$W_{\ell n}^{\tilde{\mathbf{n}}}(\tilde{\mathbf{J}})=\frac{1}{\pi }\int \mathrm{d}{\theta }_1{U}_n^{\ell }(r({\theta }_1))\cos [n_1{\theta }_1+n_2({\theta }_2-\psi )],$$

$${\theta }_1$$

$${\theta }_2-\psi$$

$${\hat{M}}_{pq}^0(\omega )$$

Full text

1 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), Gaia DR2 (Bianchini et al., 2018; Sollima et al., 2019), and spectroscopic surveys, such as the MIKiS survey (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/\sigma}$ 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. 2017, 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., 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 $N$-body simulations. Section 4 extends this analysis to the whole two-dimensional 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 distribution function 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 characterised by a mean radial potential ${\psi(r)}$, with ${(r,\vartheta,\phi)}$ the spherical coordinates. Throughout this study, we take the mean potential to be a Plummer potential (Plummer, 1911),

[TABLE]

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

[TABLE]

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\leq i\leq\pi}$. According to the Jeans theorem (Jeans, 1915), the mean-field distribution of the system is fully characterised by its DF, ${F(\mathbf{J})}$, normalised as ${\!\int\!\mathrm{d}\bm{\theta}\mathrm{d}\mathbf{J}\,F=M_{\mathrm{tot}}}$. Because the system is spherically symmetric, the orbital frequencies, ${\mathbf{\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 ${\widetilde{\mathbf{J}}=(J_{r},L)}$, and the associated reduced frequencies ${\widetilde{\mathbf{\Omega}}=\partial H(\widetilde{\mathbf{J}})/\partial\widetilde{\mathbf{J}}}$. Note that spherical systems with a spherically symmetric velocity distribution satisfy ${F=F(\widetilde{\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)

[TABLE]

with $q$ controlling the flavour and degree of velocity anisotropy of the system, and $\Gamma$ the gamma function. Here, $\mathcal{F}$ is expressed in terms111This formulation slightly simplifies Dejonghe’s by suppressing vanishing ${\Gamma(q/2)}$ terms, which are indeterminate when $q$ is an even, negative integer. of the hypergeometric function $\!\ {}_{2}F_{1}$

[TABLE]

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 the advantages of this DF is the simple forms taken by its velocity dispersions in all directions. Indeed, one has

[TABLE]

leading to the anisotropy parameter:

[TABLE]

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 (1962), a simple way of introducing rotation in a spherically symmetric system is to modify its DF as222We note that any odd fonction of $L_{z}$ could be used to generalise our choice of ${\mathrm{Sign}(L_{z})}$ to introduce rotation.

[TABLE]

with the sign function ${\mathrm{Sign}(x)}$ equal to ${-1}$ for ${x<0}$, then [math] for ${x=0}$ and $1$ for ${x>0}$. In equation (7), ${|\alpha|\leq 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,\widetilde{\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}\!\!\mathrm{d}L_{z}\,\mathrm{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, as illustrated in Figure 2.

For each orbit (i.e. each $\widetilde{\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 characterisation of the properties of these configurations will be provided in Breen et al. (2019a), here we just wish to highlight a few selected features333An implementation of the construction of such anisotropic, differentially rotating spherical models is already freely available at https://github.com/pgbreen/PlummerPlus.. In Figure 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

[TABLE]

The bottom panel of Figure 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\neq 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 Figure 3, we illustrate the 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\neq 0}$. It is possible that this discontinuity plays a role in driving the instabilities, as discussed further in Section 5.2.

3 Characterising 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, 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 ${\widehat{\mathbf{M}}(\omega)}$ (Binney & Tremaine, 2008), which in the case of ${3D}$ spherical systems is given by

[TABLE]

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

[TABLE]

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

[TABLE]

where ${\widehat{\mathbf{M}}_{0}(\omega)}$ corresponds to the system’s response in the absence of rotation, while ${\alpha\,\widehat{\mathbf{M}}_{1}(\omega)}$ captures the effects associated with rotation. Equation (11) is the master equation characterising its linear stability444We 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 ${\widehat{\mathbf{M}}_{0}(\omega)}$, the shift in eigenvalues and eigenvectors associated with the introduction of ${\alpha\ll 1}$ can be computed perturbatively.. Because of its reduced number of symmetries, we also note that the matrix ${\widehat{\mathbf{M}}_{1}(\omega)}$ is less sparse than the zero-rotation matrix ${\widehat{\mathbf{M}}_{0}(\omega)}$, and hence the numerical task is more cumbersome. In Appendix B, we detail our numerical implementation of the computation of ${\widehat{\mathbf{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 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.

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}-\widehat{\mathbf{M}}(\omega_{0}+\mathrm{i}\eta)]}$. If ever this contour encloses the origin, then the cluster supports an unstable mode with a growth rate larger than $\eta$. Figure 4 illustrates the Nyquist contours associated with our fiducial model.

The set of numerical parameters of this computation is detailed in Appendix C.1, 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 $\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 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

[TABLE]

with ${X_{p}}$ 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.

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 characterise 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

[TABLE]

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

[TABLE]

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 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 ${|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 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,\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 ${{\rm{Arg}}(A_{m}(t))}$ (​​​​$\mod{2\pi}$ to account for complete revolutions) over the same time interval. Figure 7 illustrates such measurements of $\eta$ and $\omega_{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 ${\rho^{(p)}(r,\vartheta,\phi)}$ 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 C.2, 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 ${\sim 15\%}$ discrepancy in the measured value of $\eta$, and a ${\sim 3\%}$ discrepancy in $\omega_{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.

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 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 ${\eta(\alpha,q)}$ panel (${q\!<\!0,\alpha\!\gtrsim\!\tfrac{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 ${\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 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 $\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)555Indeed, 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)..

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-biased (${q\!<\!0}$) and radially-biased systems (${q\!>\!0}$).

4.2 Instability of tangentially-biased systems: ${q<0}$

In Figure 9, 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)\mapsto\omega_{0},\eta}$ can be fitted globally by the simple polynomial expressions

[TABLE]

whose cuts are also represented in Figure 9.

In Figure 10, we represent the ${\eta(\alpha,q)}$ surface for tangentially-biased systems.

We emphasise 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 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 ${\eta(\alpha,q)}$ and ${\omega(\alpha,q)}$ surfaces in the ${q\!>\!0}$ region.

Similarly, in Figure 12, 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\!\tfrac{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 (${\alpha\!\gtrsim\!\tfrac{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 ${(\alpha,q)\mapsto\omega_{0},\eta}$ can be fitted globally by simple polynomial functions as

[TABLE]

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 ${(\alpha,q>0)}$-space.

This figure clearly highlights the transition between the two regions of instability, Fast-ROI vs. 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 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.

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_{z}$ respectively, using equations (6) and (7). They read

[TABLE]

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 reparametrisation, we present in Figure 14 a remapping of the stability manifolds in the ${(\langle\beta_{\theta}\rangle,\langle L_{z}\rangle)}$-space.

As expected the qualitative differences between Figures 8 and 14 are minor. The following simple fits represent the three manifolds

[TABLE]

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 generalised by requesting that $\eta$ is above a given threshold $\eta_{\mathrm{crit}}$. For instance, for the Slow-ROI instability, it reads ${\langle L_{z}\rangle/7+\langle\beta_{\theta}\rangle/10\geq\eta_{\mathrm{crit}}+\tfrac{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 the Lynden-Bell trick from equation (7) which may partially impact artificially the distribution of orbits near ${L_{z}=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\!\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, $Q$ was set to zero in equation (37), effectively eliminating the contribution of the discontinuity of $\mathrm{Sign}(L_{z})$ 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 $L_{z}$ scales like ${L_{z}/L}$ was also investigated. Once again, qualitatively, the three instability islands remain.

5.3 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 inertia666We note that the role played by the moment of inertia is implicit in equation (9) via the integral over $\mathbf{J}$ which can be rewritten as an integral over, say, $\Omega_{\rm ILR}$ via the corresponding moment of inertia. (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_{\mathrm{f}}=J_{r}\!+\!\tfrac{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, Figure 15 shows that most orbits have ${\zeta_{\mathrm{ILR}}\!>\!0}$, with ${\zeta_{\mathrm{ILR}}\!=\!({\partial\Omega_{\mathrm{ILR}}}/{\partial L})_{\!J_{\mathrm{f}}}}$ the cooperation parameter at the ILR, and the tumbling frequency ${\Omega_{\mathrm{ILR}}\!=\!\Omega_{\phi}\!-\!\tfrac{1}{2}\Omega_{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 $\Omega_{\mathrm{ILR}}$ is represented for various values of ${q\geq 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}}$777In such a system, all orbits precess in the same direction as they revolve., mapping them into DF$({\Omega_{\mathrm{ILR}}\mathrm{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 ${\mathrm{PDF}(\Omega_{\mathrm{ILR}})}$ will be multiplied by ${1\!-\!\alpha}$ and the right bump by ${1\!+\!\alpha}$. This rotation-induced tilting decreases the spread in ${\Omega_{\mathrm{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. 11.

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 destabilising power of this process. When the pattern speed reaches the edge of the ${\Omega_{\mathrm{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), ${\widehat{\mathbf{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 ${\langle\beta_{\theta}\rangle,\langle L_{z}\rangle}$, 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 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 ${L_{z}\!=\!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 discontinuity 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; 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 core, but it does not capture 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% (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.

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 was 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 Sklodowska-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, co-funded by CEA and CNES. SR and ALV are grateful to the School of Mathematics and Institute for Astronomy in Edinburgh as well as the Institute 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.

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)}(x)}$ and ${\rho^{(p)}(x)}$, defined as

[TABLE]

with ${U(\mathbf{x},\mathbf{x}^{\prime})=-G/|\mathbf{x}-\mathbf{x}^{\prime}|}$ 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

[TABLE]

where a given basis element is characterised by three indices, ${\ell\geq 0}$, ${|m|\leq\ell}$, and ${n\geq 0}$, and we follow the convention ${\!\int\!\mathrm{d}\vartheta\mathrm{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 w.r.t. the angle $\bm{\theta}$ as

[TABLE]

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

[TABLE]

with the reduced resonance vector ${\widetilde{\mathbf{n}}\!=\!(n_{1},n_{2})}$. Here, we introduced the rotation matrix ${R_{nm}^{\ell}(i)}$ as

[TABLE]

where the sum over $t$ is restricted to the values such that the arguments of the factorials are positive, i.e. ${t_{\rm min}\!\leq\!t\!\leq\!t_{\rm max}}$, with ${t_{\rm min}\!=\!\text{Max}[0,n-m]}$ and ${t_{\rm max}\!=\!\text{Min}[\ell-m,\ell+n]}$. In equation (22), we also introduced the Fourier-transformed “in-plane” radial coefficients ${W^{\widetilde{\mathbf{n}}}_{\ell n}(\widetilde{\mathbf{J}})}$ defined as

[TABLE]

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 $\widetilde{\mathbf{J}}$. They are given by

[TABLE]

where $\mathcal{C}$ is the integration contour going from the pericentre $r_{\mathrm{p}}$ up to the current position ${r=r(\theta_{1})}$ along the radial oscillation. Computing the radial coefficients from equation (24) 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^{[\ell^{p},-m^{p},n^{p}]}_{\mathbf{n}}=(-1)^{m^{p}}(\psi^{[\ell^{p},m^{p},n^{p}]}_{-\mathbf{n}})^{*}}$.

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. $\alpha$, the matrices ${\widehat{\mathbf{M}}_{0}(\omega)}$ and ${\widehat{\mathbf{M}}_{1}(\omega)}$ are defined as

[TABLE]

These coefficients satisfy the symmetry relations

[TABLE]

In particular, for ${\alpha=0}$, we recover the invariance of the Nyquist contours w.r.t. complex conjugation. Indeed, an immediate consequence of equation (27) is

[TABLE]

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 17, in the mildly unstable ${q\!=\!2,\alpha\!=\!0}$ case.

Owing to the Kronecker symbol, $\delta_{m^{p}}^{n_{3}}$, in equation (22), the sums over $n_{3}$ in equation (26) can be straightforwardly carried out. In equation (26), all the dependences w.r.t. $L_{z}$ are in the Fourier transformed basis elements, via the rotation matrices from equation (23). We recall that they satisfy in particular the orthogonality relation

[TABLE]

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

[TABLE]

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

[TABLE]

where the prefactor $L$ originates from the volume element ${\mathrm{d}\mathbf{J}=\mathrm{d}\widetilde{\mathbf{J}}L\mathrm{d}i\sin(i)}$. Equation (30) corresponds to the generic response matrix for a spherically symmetric and non-rotating system. One can proceed similarly to compute the matrix ${\widehat{\mathbf{M}}_{1}(\omega)}$ from equation (26). To do so, the needed gradient can be computed as

[TABLE]

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

[TABLE]

and

[TABLE]

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 relations

[TABLE]

Similarly to equation (31), let us introduce

[TABLE]

Gathering all these notations, the matrix ${\widehat{\mathbf{M}}_{1}(\omega)}$ reads

[TABLE]

Figure 18 sketches $\widehat{M}_{pq}^{0}$ and $\widehat{M}_{pq}^{1}$ for the academic case of ${n_{\mathrm{max}}=5}$ and ${\ell_{\max}=4}$ (so the corresponding block is ${6\!\times\!6}$ because indexation starts at [math]).

Both matrices are made of blocks labeled by the pair ${(\ell^{p},\ell^{q})}$, while each element of a block is labeled by ${(n^{p},n^{q})}$. Both matrices are complex (the norm of each element is shown in Figure 18) and symmetric, as can be seen by swapping $p$ and $q$ in equations (30) and (37). The off-diagonal blocks in $\widehat{M}_{pq}^{1}$ reflects interaction between different harmonics $\ell$. The ${\delta_{\ell^{p}}^{\ell^{q}}}$ in equation (30) makes $\widehat{M}_{pq}^{0}$ diagonal by blocks, while the chessboard-like structure of $\widehat{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}}}$ (equations (33) and (34)), which is non-zero only when ${\ell^{p}}$ and ${\ell^{q}}$ have the same parity.

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

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

[TABLE]

where $\rm{erf}$ (resp. $\rm{erfc}$) are the (resp. complementary) error functions, $\sigma$ is an ad hoc width, $R_{\mathrm{b}}$ is the radial scale of the basis functions. We also introduced ${u_{0}\!=\!p_{\mathrm{b}}\,n_{\mathrm{max}}\,\Delta u-1}$, with ${\Delta u}$ the step size in the ${(u,v)}$-grid, $n_{\mathrm{max}}$ the maximum radial number for the basis functions and $p_{\mathrm{b}}$ a control parameter for the step size of the ${(r_{\mathrm{p}},r_{\mathrm{a}})}$-grid where it is linearly sampled. This transformation allows for $r_{\mathrm{p}}$ to be sampled logarithmically near [math] and linearly further out, while $r_{\mathrm{a}}$ is sampled logarithmically near $r_{\mathrm{p}}$ (close to circular orbits) and linearly further out. On top of having an explicit Jacobian, equation (38) offers a smooth transition between the logarithmic and linear samplings, with ${\sigma\!=\!10^{-3}}$ the typical width of the transition region. The sampling in ${(r_{\mathrm{p}},r_{\mathrm{a}})}$ becomes linear when the radial step size reaches ${R_{\mathrm{b}}/(p_{\mathrm{b}}\,n_{\mathrm{max}})}$, hence the parameter $p_{\mathrm{b}}$ represents the minimal number of points in the $r_{\mathrm{p}}$-grid per radial oscillation for the highest order radial basis function.

Appendix B Computing ${W^{\widetilde{\mathbf{n}}}_{\ell{n}}(\widetilde{\mathbf{J}})}$

In this Appendix, we detail how the radial coefficients, ${W^{\widetilde{\mathbf{n}}}_{\ell n}(\widetilde{\mathbf{J}})}$, from equation (24) 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}^{-\widetilde{\mathbf{n}}}(\widetilde{\mathbf{J}})=W^{\widetilde{\mathbf{n}}}_{\ell n}(\mathbf{J})}$, which essentially halves the number of coefficients to evaluate.

In equation (25), 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 (24) w.r.t. $r$ rather than $\theta_{1}$. As a result, we write

[TABLE]

with the Jacobian ${\mathrm{d}\theta_{1}/\mathrm{d}r\!=\!\Omega_{1}/\sqrt{2(E\!-\!\psi(r))\!-\!J_{2}^{2}/r^{2}}}$. To compute ${\theta_{1}[r]}$ and ${(\theta_{2}\!-\!\psi)[r]}$ in equation (39), one has to compute the integrals from equation (25), with the boundaries ${\!\!\int_{\mathcal{C}}\!\mathrm{d}r^{\prime}=\!\!\int_{r_{\mathrm{p}}}^{r}\!\!\mathrm{d}r^{\prime}}$.

While all these integrals are well-defined, we note that their integrands have a diverging behaviour at the boundaries ${r=r_{\mathrm{p}},r_{\mathrm{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 w.r.t. $r$, let us therefore integrate w.r.t. an (explicit) angular anomaly, that we call $u$. First, we define ${\Sigma\!=\!(r_{\mathrm{p}}\!+\!r_{\mathrm{a}})/2}$ and ${\Delta\!=\!(r_{\mathrm{a}}\!-\!r_{\mathrm{p}})/2}$. A natural choice would be to consider the change of coordinates ${r(u)=\Sigma+\Delta\sin(\tfrac{\pi}{2}u)}$, with ${u\in[-1;1]}$, that is explicit, monotonic and satisfies ${\mathrm{d}r/\mathrm{d}u=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

[TABLE]

which is the appropriate third degree polynomial satisfying the four constraints ${f(\pm 1)=\pm 1}$ and ${f^{\prime}(\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 (39) as

[TABLE]

with

[TABLE]

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

[TABLE]

where $T_{1}$ (resp. $T_{2}$) stands for ${\theta_{1}}$ (resp. ${(\theta_{2}-\psi)}$), and are given by expressions of the form

[TABLE]

Because the integration boundaries of ${T_{1}(u)}$ and ${T_{2}(u)}$ match the ones for ${W(u)}$, one can reinterpret equation (43) 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

[TABLE]

with the initial conditions ${(W(-1),T_{1}(-1),T_{2}(-1))=(0,0,0)}$. Introducing the state vector ${\mathbf{S}=(W,T_{1},T_{2})}$, and the velocity vector ${\mathbf{s}=(w,t_{1},t_{2})}$, one can finally rewrite equation (45) under the compact form

[TABLE]

As a conclusion, under this form, computing the coefficients ${W^{\widetilde{\mathbf{n}}}_{\ell n}(\widetilde{\mathbf{J}})}$ amounts then to computing the value of ${\mathbf{S}(1)}$, by integrating forward the differential equation (46). 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 (46) 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 characterise unstable modes. Section C.1 focuses on the matrix method, while Section C.2 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)}$ .

C.1 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 (26), which is turned into a ${(r_{\mathrm{p}},r_{\mathrm{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_{\mathrm{p}},r_{\mathrm{a}})}$-grid are $r_{\min}$, $r_{\max}$, $\Delta u$, $p_{\mathrm{b}}$, $R_{\mathrm{b}}$ and were introduced in Appendix A. The parameters controlling the level of projection are $n_{1}^{\mathrm{max}}$, the maximum wave number in the Fourier decomposition over the radial angle, $n_{\mathrm{max}}$ the maximum order of radial basis functions, which also sets the size of a block in the response matrix, and $\ell_{\mathrm{max}}$ the maximum harmonic number $\ell$, which also defines the number of blocks in the response matrix (see Figure 18). We note that these two sets of parameters are not exclusive. Here, $R_{\mathrm{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_{\mathrm{max}}$ also controls the grid refinement, as can be seen in equation (38). In the following study, we made sure that varying $n_{\mathrm{max}}$ would not change the grid sampling but only the level of projection onto the radial basis.

Table 1 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 19 shows typical converging ${\eta(\alpha)}$ curves when we increase the level of precision of two of the parameters, here $p_{\mathrm{b}}$ and $n_{1}^{\max}$.

In closing, our choice of control parameters ensures that for ${\eta\!>\!0.01}$ the stability manifolds are properly mapped.

C.2 Convergence study of the $N$-body simulations

The main control parameters for the force calculation in the gyrfalcON code are $\varepsilon$ (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 $\varepsilon$ 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 20, we explore the effect of varying $\varepsilon$ (top panel) and $N$ (bottom panel).

The default value of ${\varepsilon=0.05}$ was found to give good agreement with results using lower values of $\varepsilon$, and with a direct $N$-body simulation using NBODY6 (Nitadori & Aarseth, 2012), which in effect has no softening, i.e. ${\varepsilon=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.