The Kinematic Richness of Star Clusters - II. Stability of Spherical Anisotropic Models with Rotation
Philip G. Breen, Simon Rozier, Douglas C. Heggie, Anna Lisa Varri
MMNRAS , 1–17 (2020) Preprint 9 February 2021 Compiled using MNRAS L A TEX style file v3.0
The Kinematic Richness of Star Clusters - II. Stability of SphericalAnisotropic Models with Rotation
Philip G. Breen, ★ Simon Rozier, Douglas C. Heggie, Anna Lisa Varri , School of Mathematics and Maxwell Institute for Mathematical Sciences, University of Edinburgh, Kings Buildings, Edinburgh EH9 3FD Université de Strasbourg, CNRS UMR 7550, Observatoire astronomique de Strasbourg, 11 rue de l’Université, 67000 Strasbourg, France Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
Accepted 2021 February 4. Received 2021 February 2; in original form 2020 June 29
ABSTRACT
We study the bar instability in collisionless, rotating, anisotropic, stellar systems, using 𝑁 -body simulations and also the matrixtechnique for calculation of modes with the perturbed collisionless Boltzmann equation. These methods are applied to sphericalsystems with an initial Plummer density distribution, but modified kinematically in two ways: the velocity distribution istangentially anisotropic, using results of Dejonghe, and the system is set in rotation by reversing the velocities of a fraction ofstars in various regions of phase space, à la Lynden-Bell. The aim of the 𝑁 -body simulations is first to survey the parameter space,and, using those results, to identify regions of phase space (by radius and orbital inclination) which have the most importantinfluence on the bar instability. The matrix method is then used to identify the resonant interactions in the system which havethe greatest effect on the growth rate of a bar. Complementary series of 𝑁 -body simulations examine these processes in relationto the evolving frequency distribution and the pattern speed. Finally, the results are synthesised with an existing theoreticalframework, and used to consider the old question of constructing a stability criterion. Key words: galaxies: kinematics and dynamics – galaxies: star clusters: general
The analysis of the stability properties of self-gravitating rotatingequilibria is a classical problem in fluid and stellar dynamics (Chan-drasekhar 1969), starting with the study of the ellipsoidal figures ofequilibrium. The sequence of the uniformly rotating axisymmetricMaclaurin spheroids, as parameterised by the ratio of the rotationalkinetic energy to the gravitational energy 𝑇 /| 𝑊 | , exhibits a bifurca-tion point (at 𝑇 /| 𝑊 | = . 𝑇 /| 𝑊 | which are greater than such a threshold are secularly unstable with respect to global bar modes. Further along the Maclau-rin sequence (when 𝑇 /| 𝑊 | > . dynamical instability withrespect to a bar mode sets in.In the theory of rotating stars (e.g., see Tassoul 1978 for a com-prehensive introduction), the investigation of the properties of self-gravitating rotating fluid bodies has been generalized to the case ofconfigurations with nonuniform density, especially polytropic flu-ids with solid-body (e.g., see James 1964; Stoeckly 1965; Lai et al.1993) and with differential rotation (e.g., see Ostriker & Mark 1968;Tohline et al. 1985; Hachisu 1986; Pickett et al. 1996; New et al.2000).On the other hand, in the study of rotating stellar dynamical sys-tems, stability analyses have been limited mainly to the context ofdisks. Hohl (1971) and Kalnajs (1972) provided the numerical andanalytical evidence that uniformly rotating disks can be strongly un-stable with respect to global bar modes, if the ordered kinetic energy ★ E-mail: [email protected] dominates the energy budget of the system. Ostriker & Peebles (1973)extended this numerical investigation to the case of differentially ro-tating disks and conjectured the condition 𝑇 /| 𝑊 | < . ± .
02 as an empirical necessary (but not sufficient) criterion for the dynamicalstability of any rotating stellar system with respect to global bar-likemodes.Over the years, a number of studies, mostly based on N-bodysimulations, have been carried out to investigate the stability of stellardisks (e.g., see Hohl 1976; Sellwood 1981), so that the nature ofthe rotational instability with respect to bar-like modes was indeedconfirmed to be dynamical. Few attempts at providing a physicalinterpretation of the Ostriker & Peebles (1973) criterion have beenmade (Vandervoort 1982); several counterexamples or alternativestability criteria have also been suggested (e.g., see Zang & Hohl1978; Berman & Mark 1979; Aoki, Noguchi & Iye 1979; Efstathiouet al. 1982; Evans & Read 1998; Athanassoula 2008).Yet, the stability properties of differentially rotating spheroidal stellar systems have been rarely explored (see Palmer, Papaloizou &Allen 1990; Kuijken & Dubinski 1994; Sellwood & Valluri 1997) andthe connection with the corresponding fluid systems is only partiallyunderstood (see Christodoulou et al. 1995). Some exceptions arerepresented by the collisionless counterparts of the uniformly rotatingpolytropes (Vandervoort 1980), Riemann ellipsoids (Vandervoort &Welty 1982), and Maclaurin spheroids (Vandervoort 1991), for whichanalytical studies have demonstrated that the dynamical instabilitywith respect to bar mode sets in at lower values of rotation withrespect to the Maclaurin fluid sequence (e.g., at 𝑇 /| 𝑊 | = . © a r X i v : . [ a s t r o - ph . GA ] F e b Breen, Rozier, Heggie & Varri. of the stability of differentially rotating fluids, kindled by the sur-prising discovery by Centrella et al. (2001) of an unstable 𝑚 = 𝑛 = .
33, and 𝑇 /| 𝑊 | ≈ .
14. Since then, numerical studies haveconfirmed that 𝑚 = , 𝑇 /| 𝑊 | as lowas 0 .
01 (Shibata et al. 2002; Saijo et al. 2003; Karino & Eriguchi2003). The study of the stability of differentially rotating sphericalshells (Watts et al. 2003, 2004) suggests that the unstable modeswithin a low- 𝑇 /| 𝑊 | differentially rotating configuration are charac-terised by corotation within the system. In particular, Watts et al.(2005) argued that differential rotation must be sufficiently strong inorder for instability to occur. Numerical investigations of the stabil-ity of differentially rotating polytropes by Saijo & Yoshida (2006)and Ou & Tohline (2006) strengthened this picture, which has nowbeen formalised also in an analytical setting (Saijo & Yoshida 2016;Yoshida & Saijo 2017).In a broader stellar dynamical context, it has been known for overa century that corotation plays a central role in the stability of shearflows (e.g., see Orr 1907; Lin 1955); in the case of self-gravitatingstellar disks, this intuition has offered a key guiding principle sincethe beginning of the development of the density wave theory ofspiral structure (e.g., see Binney & Tremaine 2008, Bertin & Lin1996 and references therein), especially in the context of variouspossible mechanisms underpinning the formation of a central bar indisk galaxies (e.g., see pioneering studies by Lynden-Bell & Kalnajs1972; Mark 1977; Toomre 1981).Leaving aside for a moment the role of the angular momentum,it is worth remembering here that a number of instructive investiga-tions have been focussed on the behaviour of spherical , non-rotatingstellar systems in a purely collisionless setting. Also in this case,the analogy with simple fluid systems has offered a fruitful path-way: a numerical stability analysis of the collisionless counterpartsof classic polytropes has been conducted by Hénon (1973), whilean analytic criterion for the stability of the case of polytropic index 𝑛 = / MNRAS , 1–17 (2020) tability of Rotating Anisotropic Models of angular momentum in collisional stellar dynamics (Part III, Breenet al., in preparation).The present article is organised as follows. Section 2 describesthe classes of equilibria adopted as initial conditions of our N-bodysimulations and outlines the methodology employed for the studyand characterisation of the relevant unstable modes. The results ofour stability analysis are presented in Section 3, together with a studyof the spatial-temporal structure of the instability. Next, in Section 4,we investigate the role of resonances, supplementing the analysis ofour N-body models with a comparative study conducted with the aidof the matrix method (as customised to rotating spheres by Rozier, etal. 2019). Finally, in Section 5 we discuss our results within existingtheoretical mechanisms for bar formation and the broader context ofglobal stability criteria for rotating, anisotropic stellar systems. 𝑁 -body methods For the study of collisionless stability in stellar systems, the two mostimportant methods are 𝑁 -body techniques and matrix methods. Thelatter deal only with linear stability, and tend to be less flexible, butprovide a direct insight into the dynamics, being naturally expressiblein terms of action-angle variables (or other invariants). We, therefore,adopt this approach in Sec.4.2, where we aim to study the role ofvarious resonances in the observed instabilities. In order to explorea wide range of initial conditions (especially in the rotation profile),however, the results we present first (Sections 3.1–3.3) adopt parti-cle methods. The following paragraphs consider these two kinds oftechniques in turn.Generally, we avoid direct-summation 𝑁 -body methods, becausethe need to complete many simulations in a reasonable time tendsto restrict systems to modest 𝑁 , when stability may be modifiedby collisional processes. It has even been said that they “shouldbe regarded as methods of last resort” (Sellwood 1997). Thereforewe mainly adopt an efficient (hierarchical) method which makes itfeasible to carry out simulations with larger particle numbers, andemploys softening to further reduce collisionality; for these reasonsour method of first resort in this work is gyrfalcON (Dehnen 2002).This is used to generate the main bulk of the numerical results, whichare in the earlier parts of Sec.3. Elsewhere, however, in Sections 3.4and 4.3, we do resort to use of the collisional code NBODY6 (Nitadori& Aarseth 2012).The matrix method dates back to the work of Kalnajs 1977 (thoughhis emphasis was on disks) and to the study of non-rotating sphericalsystems by Polyachenko & Shukhman (1981). For such systems theformulation we employ (Rozier, et al. 2019, see, for example, Eq.(4)in the present paper) would be identical to that in the Appendix ofPolyachenko & Shukhman (1981), except for notation, which in ourcase is derived from Hamilton et al. (2018). Our choice of basisfunctions (of biorthogonal density-potential pairs) can be found inHernquist & Ostriker (1992) and relies on ultraspherical polynomials.In the case of rotating systems our formulation of the matrix equationmay also differ a little from those found in the literature (e.g. Palmer1994, eq.(8.18)), not only because of notation, but also through thechoice of distribution function, though such differences are minor.Throughout the paper (even in sections not devoted to 𝑁 -bodysimulations), units are such that 𝐺 , the total mass and virial radiusare unity (i.e. Hénon Units). The initial conditions we adopt all have the same density profile, thatof a Plummer model, and the velocity distribution is constructed ina series of steps:(i) We begin with a model from the series constructed by De-jonghe (1987). These are parametrised by a single dimensionlessparameter, 𝑞 , in the range (−∞ < 𝑞 ≤ ) . Increasing positive val-ues correspond to models with increasing radial anisotropy, whileincreasingly negative values correspond to those with increasinglytangential anisotropy. The Einstein-Plummer model, with all stars oncircular orbits, corresponds to the limit 𝑞 → −∞ .(ii) The velocities are adjusted so that a fraction 𝛼 of the particleswith negative azimuthal velocity 𝑣 𝜙 (in spherical coordinates 𝑟, 𝜃, 𝜙 )are changed from negative to positive in a certain region of phasespace, as described below. Thus 𝛼 = 𝛼 = 𝑣 𝜙 > <
0, respectively. (Thedemon can be applied in reverse, as in (b) below.)(a)
Basic models:
In the simplest series of models, the demon isallowed to work only on a certain fraction of the mass (0 ≤ 𝛾 ≤ 𝛾 = . 𝛾 = 0 gives a Plummer-Dejonghe model, while when 𝛾 = 𝛼 =
1, all particles are given prograde rotation. These models areintended to give insight into the range of radii which are active inthe instabilities observed in them.(b)
High-shear models:
In the next series of models, the aboveprocedure is applied to all particles (or rather to a fraction 𝛼 ofthem) in the most bound fraction 𝛾 , but in the rest, i.e. the leastbound fraction 1 − 𝛾 of them, the Lynden-Bell demon acts inreverse, reversing 𝑣 𝜙 for all prograde stars (or a fraction 𝛼 ofthem). These models are intended to maximise the shear in therotation profile (Sec.2.3), as this may be a factor in the existenceof instabilities.(c) Low-inclination models:
Third, we consider a series of mod-els in which the demon works only on stars of low inclination 𝑖 < 𝛾 𝜋 /
2, where 𝑖 is defined in the sense of the inclination of theorbital plane to the equatorial plane 𝑧 =
0, so that 0 ≤ 𝑖 ≤ 𝜋 /
2. Amotivation here is to investigate to what extent an instability canbe thought of as being driven by motions in a thick conical disk.(d)
High-inclination models:
Lastly we have a set of models inwhich the Lynden-Bell demon works only on stars of high incli-nation 𝑖 > 𝛾 𝜋 /
2. Such models may give additional indicationsof the geometrical and kinematic characteristics of orbits whichcontribute strongly to instability.These are the only initial models considered in this paper, butmany variants are possible. For example, in the high-shear models,the axis of rotation of the least bound stars may be chosen differentlyfrom that of the most bound stars. Many such models can be exploredby setting up initial conditions with the software
PlummerPlus . There are several ways in which the rotation of these models canbe measured. The first is a rotation profile, which can be defined in https://github.com/pgbreen/PlummerPlus MNRAS , 1–17 (2020) Breen, Rozier, Heggie & Varri.
Figure 1.
Rotation curves of six “basic" models with 𝑞 = − , 𝛼 = 𝛾 = ( . )
1. Here 𝛾 is the fraction of stars, ordered by energy,below which the Lynden-Bell demon is applied. The ordinate is the meanazimuthal velocity in a cylindrical shell of radius 𝑅 (the abscissa). Each of the200 shells contains 500 particles. The six plots correspond to application ofthe Lynden-Bell demon (with various values of 𝛾 ) to the same realisation ofa Plummer-Dejonghe model; therefore the fluctuations in successive modelsare correlated. different ways, e.g. in spherical or cylindrical geometry, of whichwe concentrate on the latter. It is easily calculated from any modelby first calculating the average rotation speed ¯ 𝑣 𝜙 ( 𝑅 ) = (cid:104) 𝑣 𝜙 (cid:105) asthe average azimuthal velocity in a cylindrical shell of cylindricalradius 𝑅 , and a thickness chosen as a compromise between sufficientspatial resolution and sufficient statistical reliability. From an 𝑁 -bodysimulation the corresponding quantity is easy to compute on a gridof 𝑅 . Several examples are given below (Sec.2.4).For a global measure of rotation we compute the rotational kineticenergy, but for this purpose define the average rotation speed on acircular ring of given 𝑅 and 𝑧 and, again, suitably chosen bin size.Then the rotational kinetic energy is defined to be 𝑇 = ∫ 𝜌 ( 𝑅, 𝑧 ) ¯ 𝑣 𝜙 ( 𝑅, 𝑧 ) 𝜋𝑅𝑑𝑅𝑑𝑧, (1)the integral being evaluated in practice by summing over the chosengrid. For our purpose a fixed grid in 𝑅 with 100 bins, and a gridin 𝑧 adapted to contain approximately equal numbers of particles,was adopted. Of particular importance in this field, because of thepioneering work of Ostriker & Peebles (1973), is the dimensionlessparameter 𝑇 / 𝑊 , where 𝑊 is the binding energy of the system. Fig.1 illustrates various rotation curves for the set of “basic" models(Sec.2.2). In this figure the uppermost (blue) curve shows the rotationcurve of a model where the Lynden-Bell demon is applied to allstars, whereas the lowest shows the rotation curve of the originalPlummer model (anisotropy in all the models being modified bysetting 𝑞 = − 𝑁 -body model. The green curve, correspondingto 𝛾 = .
2, shows the result when the demon is applied only tothe 20% most bound particles. Other choices of 𝛼 simply alter theordinate in proportion to 𝛼 .Fig.2 gives similar plots for the “high-shear" (counter-rotating) Figure 2.
Rotation curves of six “high-shear" models with 𝑞 = − , 𝛼 = 𝛾 = ( . )
1. As Fig.1, except that here 𝛾 is the fraction ofstars, ordered by energy, below which the Lynden-Bell demon is applied inthe sense of prograde rotation, while above that energy it is applied in such away as to give retrograde rotation. models, in which 𝛾 is the fraction of stars, ordered by energy, belowwhich the Lynden-Bell demon is applied to give prograde rotation,while it is applied to the less bound stars to give counter-rotation.Qualitatively the rotation curves look very similar, but the differencesbetween the maximum and minimum (largest retrograde) velocitiesare about twice as large.Next, we can consider the effect of the inclination of the rotatingorbits, by applying the Lynden-Bell demon only to orbits of lowinclination; specifically, to stars such that 𝑖 < 𝛾 𝜋 /
2. The resultingrotation curves are illustrated in Fig.3. At small cylindrical radius thisprocedure greatly reduces the amount of rotation (unless 𝛾 is verylarge), since most stars at small cylindrical radius must lie on high-inclination orbits. The effect at large radius is much smaller, unlessthe critical value 𝛾 is very small. In a complementary fashion wecan also choose to apply the demon only to high-inclination orbits,such that 𝑖 > 𝛾 𝜋 /
2, which leads to the rotation curves in Fig.4.
Our method of searching for a bar, in a spherical model which isrotating about the 𝑧 − axis, is described in Rozier, et al. (2019), alongwith studies of its accuracy . Briefly, a bar mode in the 𝑥, 𝑦 plane ismeasured by computing the following integral, which uses cylindricalpolar coordinates: 𝐶 𝑚 ( 𝑡 ) = ∫ ∞ 𝑅 d 𝑅 ∫ ∞−∞ d 𝑧 ∫ 𝜋 d 𝜙 𝜋 𝜌 ( 𝑅, 𝑧, 𝜙, 𝑡 ) e − i 𝑚𝜙 , (2)where 𝑚 = to | 𝐶 𝑚 ( 𝑡 )| over a suitable time interval. Actu-ally this is a relatively insensitive measure for non-rotating systems( 𝛼 = 𝑥, 𝑦 plane. Therefore non-rotating systems ( 𝛼 =
0) can be mea-sured by computing the corresponding integrals involving spherical See also the caption to Fig.5. A logistic fit is used in Appendix A.MNRAS , 1–17 (2020) tability of Rotating Anisotropic Models Figure 3.
Rotation curves of six “low-inclination" models with 𝑞 = − , 𝛼 = 𝛾 = ( . )
1. As described in the text, the Lynden-Bell demonis applied only to low-inclination orbits, with 𝑖 < 𝛾 𝜋 /
2. The ordinate is themean azimuthal velocity in a cylindrical shell of radius 𝑅 (the abscissa). Eachof the 200 shells contains 500 particles. Figure 4.
Rotation curves of six “high-inclination" models with 𝑞 = − , 𝛼 = 𝛾 = ( . )
1. The Lynden-Bell demon is applied only tohigh-inclination orbits, with 𝑖 > 𝛾 𝜋 /
2. The ordinate is the mean azimuthalvelocity in a cylindrical shell of radius 𝑅 (the abscissa). Each shell contains500 particles. harmonics of order 𝑙 =
2, but this is not what is done in the presentstudy.
We begin (Sec.2.2) with the models in which prograde rotation isadded at energies less than (more bound than) a certain fraction 𝛾 of mass. Let us consider first the fiducial case 𝑞 = −
6. Fig.5 shows thegrowth rate 𝜂 as a function of the Ostriker-Peebles parameter 𝑇 / 𝑊 ,where 𝑇 is the rotational (mean flow) part of the kinetic energy. Notethat results for 𝛾 = Figure 5.
Growth rate versus Ostriker-Peebles parameter for six sets of modelswith 𝑞 = − , ≤ 𝛼 ≤ 𝛾 . Note that two realisations have beencomputed for the case 𝛾 =
1, which gives an impression of the accuracywith which the growth rate has been estimated. the sense that 𝜂 is not significantly non-zero), according to Rozier,et al. (2019). Results for 𝛾 = 𝛼 , to give an indication of thenumerical uncertainty in the measured growth rate. Each sequenceof symbols corresponds to a given value of 𝛾 , and one moves to theright within a given sequence as 𝛼 increases. This increase in 𝜂 issatisfactorily fitted by a cubic power law in 𝑇 / 𝑊 , or equivalently 𝜂 ( 𝑞, 𝛾 , 𝛼 ) = 𝜂 𝑚𝑎𝑥 ( 𝑞, 𝛾 ) 𝛼 , (3)a remarkably steep dependence.Because of the scaling with 𝑇 / 𝑊 , we shall now only consider thefunction 𝜂 𝑚𝑎𝑥 ( 𝑞, 𝛾 ) , results for which are shown in Fig.6. Resultsof two models are given at 𝛾 =
1, again to give an impression of thereproducibility of the results for independent initial conditions. Eventhe large fluctuations (such as the maximum at 𝛾 = . 𝛾 increases. Rather, the growthrate saturates when roughly the inner half of the system (ordered byenergy) has been set to rotate in one sense . Turning back to Fig.1,in particular the rotation curves for 𝛾 = . , . 𝑟 ∼ .
5, make little difference to the growth rate. By contrast therotation curves for 𝛾 = . . 𝑟 ∼ .
5, as do the growth rates (Fig.6). This region ( 𝑟 < ∼ . 𝑣 𝜙 ( . ) . This plot shows that this value of ¯ 𝑣 𝜙 is a roughly linearpredictor of the growth rate, though for very negative values of 𝑞 itdoes not predict well the saturation of the growth rate at the largest Actually there is some evidence from this figure that the saturation levelis a little lower than the peak level, as one notices particularly for the mostnegative values of 𝑞 . MNRAS000
5, as do the growth rates (Fig.6). This region ( 𝑟 < ∼ . 𝑣 𝜙 ( . ) . This plot shows that this value of ¯ 𝑣 𝜙 is a roughly linearpredictor of the growth rate, though for very negative values of 𝑞 itdoes not predict well the saturation of the growth rate at the largest Actually there is some evidence from this figure that the saturation levelis a little lower than the peak level, as one notices particularly for the mostnegative values of 𝑞 . MNRAS000 , 1–17 (2020) Breen, Rozier, Heggie & Varri.
Figure 6.
Growth rate for the maximally rotating model at each 𝛾 , for severalvalues of 𝑞 . We have avoided the isotropic case 𝑞 =
0, as it is known (Rozier,et al. 2019) that any instabilities here could belong to a different family fromthose at smaller (more negative) values of 𝑞 . The result of two independentsimulations is given for 𝛾 = Figure 7.
Growth rate for the maximally rotating model at each 𝛾 , for severalvalues of 𝑞 , plotted against the value of the rotation speed at radius 𝑅 = . values of ¯ 𝑣 𝜙 ( . ) . Incidentally, this value was conveniently available,but values at other neighbouring radii might well be equally useful,as illustrated just below.The maximum growth rate occurs when 𝛼 =
1, i.e. the Lynden-Bell demon acts at its maximum level. But now the entire dataset (allvalues of 𝑞 ≤ , 𝛼 and 𝛾 ; see Fig.8) can be grasped in outline bymeans of the foregoing analysis. The axes correspond to the Ostriker-Peebles stability parameter and the azimuthal streaming speed at thevirial radius, and the size of the symbol encodes the growth rateof the bar. The region of the diagram containing the fastest growthrate is particularly interesting: although the classical Ostriker-Peeblescriterion (for flattened systems) 𝑇 / 𝑊 > ∼ .
14 corresponds rather wellwith the boundary of this region when ¯ 𝑣 𝜙 ( ) is large, we now findstrongly unstable models when 𝑇 / 𝑊 is much smaller, provided that¯ 𝑣 𝜙 ( ) is small. These are models in which 𝛾 is sufficiently smallthat the azimuthal speed has dropped almost to zero by the virial Figure 8.
Scatterplot of all results on models with prograde rotation, in theplane of 𝑇 / 𝑊 , ¯ 𝑣 𝜙 ( ) , i.e. the Ostriker-Peebles parameter and the azimuthalspeed at the virial radius. The growth rate is coded in the size of the plottedpoints; the smallest symbols indicate growth rates less than 0.05. radius (cf. Fig.1). In fact the evident, roughly horizontal groupingsin this diagram correspond to the two or three highest values of 𝛾 included in our numerical survey. For this sequence of models, we assume that the scaling with 𝑇 / 𝑊 is similar to that in the purely prograde case, and proceed directlyto the results for 𝛼 =
1. The growth rate is shown in Fig.9. Qual-itatively the results behave quite similarly to those for the purelyprograde models, and even the maximum growth rates (which occurfor models with smaller, i.e. more negative, values of 𝑞 , and then atintermediate values of 𝛾 (cid:39) .
5) are not as dissimilar as the moreextreme gradient in the velocity curve (compare Figs.1, 2) mightlead one to expect. Nevertheless one notices even more strongly thatthe maximum growth rate occurs at intermediate values of 𝛾 (cf.footnote 4); furthermore the growth rate there is even larger than inFig.6. These remarks suggest that the maximum growth rate dependson the rotation of the outer (least bound) parts of the system: as 𝛾 increases beyond 0.5, these parts start to rotate, and the growth rateis depressed slightly, while as 𝛾 increases beyond 0.5 (Fig.9) theouter regions switch from retro- to prograde rotation, and the drop ingrowth rate is more pronounced and starts from a higher maximum.In Fig.9 there are significant changes in the rate of growth atthe left, around 𝛾 = .
2, where the values of 𝜂 𝑚𝑎𝑥 have a verypronounced minimum, increasing again as 𝛾 decreases further. Infact for 𝛾 = 𝛾 =
1, as these two cases have the samerotation curve except for the sense of rotation.If we contrast the great difference in 𝜂 𝑚𝑎𝑥 between 𝛾 = 𝛾 = . 𝑅 =
1, it isclear that ¯ 𝑣 𝜙 ( ) (coupled with 𝑇 / 𝑊 ) cannot be as simple a predictorof the growth rate as in the case of the prograde models (Fig.8). InFig.10 we show all results for the counter-rotating models, exceptthat now the ordinate is the absolute value of ¯ 𝑣 𝜙 ( ) . The resultis qualitatively similar to that for the prograde models, but morenoticeable here is the fact that there are some very rapidly rotatingmodels (high 𝑇 / 𝑊 and | ¯ 𝑣 𝜙 ( )| ) with relatively low growth rates. We MNRAS , 1–17 (2020) tability of Rotating Anisotropic Models Figure 9.
Growth rate for the maximally rotating (and counter-rotating) modelat each 𝛾 , for several values of 𝑞 . For the reason stated in the text, the resultsfor 𝛾 = 𝛾 = Figure 10.
Scatterplot of all results on models with pro- and retrograderotation, in the plane of 𝑇 / 𝑊 , | ¯ 𝑣 𝜙 ( ) | , i.e. the Ostriker-Peebles parameterand the absolute value of azimuthal speed at the virial radius. The growthrate is coded in the size of the plotted points; the smallest symbols indicategrowth rates less than 0.05. now investigate the rotation curves of such models to try to identifydistinctive features of their rotation curves.In fact the models with low growth rate, high 𝑇 / 𝑊 and high | ¯ 𝑣 𝜙 ( )| tend to have small (very negative) 𝑞 and 𝛾 ∼ .
2. Such models arealso very distinctive in Fig.9, as discussed above. What distinguishesthe rotation curve for these models, as can be seen in Fig.2, is thatthe interior radii of a 𝛾 = . 𝑣 𝜙 ( . ) is small.Therefore, to separate such models better than in Fig.10, we simplymultiply the abscissa by | ¯ 𝑣 𝜙 ( . )| . The result is shown in Fig.11.There is now a much better separation of slowly and rapidly growingmodes, and this choice of parameters has comparable success in thebasic models of Sec.3.1. These results confirm our finding (in that Figure 11.
As Fig.10, except that the abscissa is multiplied by the meanazimuthal speed at radius 𝑅 = .
2, to separate models which are rotatingrelatively slowly at small radii. The growth rate is coded in the size of theplotted points; the smallest symbols indicate growth rates less than 0.05. subsection) that the kinematics of the region inside 𝑅 (cid:39) . 𝛾 = 𝑞 = − Ω 𝑏 (cid:39) − . . The pattern speed reduces in magnitude to-0.5 at 𝛾 = .
1, but is about +0.7 (i.e. prograde) throughout therange 0 . ≤ 𝛾 ≤
1. At 𝛾 = . 𝑡 (cid:39)
40 aweak bar grows with a fast pattern speed around Ω 𝑏 = − .
7, but thisgives way by 𝑡 (cid:39)
60 to a bar with faster growth and Ω 𝑏 (cid:39) − . More briefly, we consider the maximal growth rate (i.e. for 𝛼 =
1) formodels in which rotation is given only to orbits of low inclination,i.e. those with rotation curves exemplified by Fig.3, or only to orbitsof high inclination (Fig.4). The results are shown in Fig.12, withthe low-inclination series on the right. As 𝛾 decreases from 1 toabout 0.8, the growth rate drops by about a factor of 2, showingthat stars with inclination above about 70 ◦ are very important inthe instability. (This does not mean that stars of low inclinationare irrelevant, but they are not sufficient of themselves for rapidlygrowing instability.) A complementary conclusion is drawn from thehigh-inclination models, whose growth rate drops by a comparablefactor if 𝛾 increases from 0 to about 0.3, i.e. if one excludes orbitswith inclination less than about 30 ◦ . Thus the presence of both high-and low-inclination stars is necessary for rapid growth, but neithergroup by itself is sufficient. Measurements reported in this paragraph were taken in a time intervalaround 7 < 𝑡 <
21, from the time of the clear emergence of the bar up to atime which is about half of the time when it reaches its maximum amplitude.The values are thought to be accurate to better than 10%.MNRAS000
21, from the time of the clear emergence of the bar up to atime which is about half of the time when it reaches its maximum amplitude.The values are thought to be accurate to better than 10%.MNRAS000 , 1–17 (2020)
Breen, Rozier, Heggie & Varri.
Figure 12.
As Fig.9, except that the Lynden-Bell demon is applied at only lowor high inclinations, as determined by the parameters 𝛾 and 𝛾 respectively. In this subsection we focus on the basic case 𝑞 = − , 𝛼 = 𝛾 = 𝑞 = −∞ ). Our aim is to present brieflysome information on the spatial extent and shape of the bar. For othervalues of 𝑞 , representations of the density are shown in Rozier, et al.(2019, figs.5,6 for 𝑞 = −
6, and fig.10 for other 𝑞 < 𝑁 -body simulations with 𝑁 = NBODY6 , the bar is detectable (by Fourier analysis) only insidethe cylinder
𝑅 < ∼ . , | 𝑧 | < ∼ .
9. Fig.13 displays the projected densityin the plane 𝑧 = 𝑡 =
20 (in this model the bar reaches itsmaximum amplitude at about 𝑡 = 𝑁 -body model at time 𝑡 =
20, the core radius has decreased from an initial value of about0.32HU to a value of about 0.26HU, presumably because of earlycore collapse. Possibly for this reason, the projected density at thecentre of the bar is above zero.
The results that follow are based on a dedicated set of direct 𝑁 -body models with 𝑁 = 𝑞 = −
16. They are basic models of the type discussed in Sec.3.1 with 𝛾 =
1, i.e. the Lynden-Bell demon is applied everywhere. Our aimis to discuss the dependence of the time scales of the bar instabilityon the rotation parameter 𝛼 , with a view to informing the discussionof Section 5.1.2. Detailed discussion of the methods of analysing theruns are given in Appendix A. Results are given in Fig.14.While details of the analysis and interpretation of this plot aregiven in later sections and the appendix, explanation is required herefor the fact that any pattern speed at all has been measured for the non-rotating case 𝛼 =
0. We emphasise that no growth could be measured
Figure 13.
Projected density, in the plane 𝑧 = 𝑁 -body modelwith 𝑁 = 𝑡 = 𝑡 =
0, and therefore shows only the perturbation. Thedata was gridded with a pixel size 0 . HU, and contours were smoothed byinterpolating to a grid 9 times finer. The scale gives the number per pixel.
Figure 14.
Growth rate and pattern speed of direct 𝑁 -body models with 𝑁 = , 𝑞 = − , 𝛾 = 𝛼 . For 𝛼 = , . 𝛼 = . 𝛼 . For 𝛼 = Ω − Ω / with confidence; all that is shown in Fig.14 is an upper limit. Butthe analysis of the 𝑚 = 𝐶 (see Eq.(2)), usingthe autocorrelation of arg ( 𝐶 ) , clearly demonstrated the presence ofoscillation with the stated pattern speed. The figure supports the viewthat the bar instability in models with 𝛼 > MNRAS , 1–17 (2020) tability of Rotating Anisotropic Models Table 1.
Resonances in the leading-order perturbations of epicyclic motionsin a rotating spherical model. Names which are not in common use are givenin slanted type. Here the orbital inclination 𝑖 is used in its conventional sense,with retrograde orbits having 𝜋 / < 𝑖 ≤ 𝜋 . Ω 𝑏 = 𝑖 -dependence Name ( 𝑛 , 𝑛 ) 𝑖 tumbling (0,0) Ω ( + cos 𝑖 ) corotation ( , ) Ω − 𝜅 / ( + cos 𝑖 ) ILR (− , ) Ω + 𝜅 / ( + cos 𝑖 ) OLR ( , )± 𝜅 / 𝑖 epicyclic (± , )− Ω ( − cos 𝑖 ) reverse corotation ( , − )− Ω + 𝜅 / ( − cos 𝑖 ) reverse ILR ( , − )− Ω − 𝜅 / ( − cos 𝑖 ) reverse OLR (− , − ) The role of resonances in the occurrence of bar instabilities in disksis sufficiently well known to be textbook material (e.g. Binney &Tremaine 2008). Less familiar, however, is their role in sphericalsystems, despite the substantial literature referred to in Sections 1 and2.1. In this introductory subsection we focus on (very) tangentiallyanisotropic systems, where the radial motions are small, and we canemploy the epicyclic approximation. Then the main resonances canbe identified by a quite elementary calculation (Appendix B); it leadsto the results given in Table 1, which we now explain.We suppose there is a bar perturbation, whose pattern rotates withangular speed Ω 𝑏 . Then the familiar corotation resonance (line 2) isgiven by Ω 𝑏 = Ω , where Ω is the circular angular velocity. Similarlythe Inner and Outer Lindblad Resonances (ILR, line 3; OLR, line 4)are defined by Ω 𝑏 = Ω ∓ 𝜅 /
2, where 𝜅 is the epicyclic frequency.What is more unfamiliar is their inclination-dependence; and whatwe have termed reverse resonances would be the normal corotationand Lindblad resonances for a bar with retrograde pattern speed (asone can see from the 𝑖 -dependence). Finally, there are what we havetermed epicyclic and tumbling resonances, which are strongest athigh inclinations (near 𝜋 / Ω 𝑏 = 𝑛 Ω + 𝑛 Ω , where 𝑛 , 𝑛 are integers, Ω isthe frequency of radial motions, and Ω is the mean longitudinalfrequency, which generalise 𝜅 and Ω , respectively. These numberpairs are included in the last column of the table.The simplest situation in which to visualise the possible location ofthese resonances is the limiting case 𝑞 → −∞ (the “Einstein sphere"),when the frequencies Ω , Ω become 𝜅, Ω . Combinations of thesefrequencies are plotted in Fig.15. Also shown is the observed patternspeed (Sec.3.4.2) found for the model with 𝑞 = −
16 in the cases 𝛼 = . ,
1. Here the anisotropy is sufficiently tangentially biased thatthe Einstein sphere is approximately representative.In a later section (Sec.4.3) we shall present some evidence of theinfluence of resonances on the behaviour of particles in 𝑁 -body mod-els. Here, for one model, we simply present in Fig.16 the distributionof particles in the space of the frequencies Ω , Ω , where their prox-imity to the various resonances can be appreciated visually. For theinitial conditions of a model with 𝑁 = 𝑞 = −
16 and 𝛼 = 𝛾 = Ω , were calculated byintegrating the equations of motion in the Plummer potential.The distribution is narrow because the anisotropy is quite extreme( 𝑞 = −
16) (in the limit 𝑞 → −∞ the distribution would have beenone-dimensional). The Inner Lindblad Resonance and its retrogradecounterpart appear to involve no particles, but we shall see (Sec.4.2.2) Figure 15.
Resonances related to circular motions in a Plummer potentialperturbed by a bar. Also shown is the approximate pattern speed of a bar in arotating model with 𝑞 = − , 𝛼 = Figure 16.
Distribution of frequencies Ω , Ω at time 15 in an 𝑁 -bodymodel with 𝑁 = , 𝑞 = − , 𝛼 = , 𝛾 =
1, along with the loci offive resonances. At this time the model exhibits a growing bar (Sec.4.3),and the pattern speed was determined in its early growth phase. OLR/ILR:Outer/Inner Lindblad Resonance, CR: Corotation Resonance, ER: epicyclicresonance (Table 1), rILR: reverse ILR. that the former still appears to play a central role in the dynamics, andpresent below a possible reason for this. There are, however, particleslying close to the other three resonances, and empirical evidencefor their influence is studied in Sec.4.3. We can see immediately,however, that these particles are confined to the low-frequency partof the distribution, and this means that they affect stars at relativelylarge radius. With the exception of the “tumbling" resonance, theremaining resonances in Tab.1 appear to play no role.In this section we identified resonances by examining the Fourierdecomposition of the perturbing potential. In first-order orbital per-turbation theory, a given Fourier term gives rise to perturbationscontaining a resonant denominator 𝜔 − 𝑛 Ω − 𝑛 Ω , where 𝜔 is thefrequency of the explicit time-dependence of the perturbing potential.If resonance is strictly defined as the vanishing of this denominator MNRAS000
1, along with the loci offive resonances. At this time the model exhibits a growing bar (Sec.4.3),and the pattern speed was determined in its early growth phase. OLR/ILR:Outer/Inner Lindblad Resonance, CR: Corotation Resonance, ER: epicyclicresonance (Table 1), rILR: reverse ILR. that the former still appears to play a central role in the dynamics, andpresent below a possible reason for this. There are, however, particleslying close to the other three resonances, and empirical evidencefor their influence is studied in Sec.4.3. We can see immediately,however, that these particles are confined to the low-frequency partof the distribution, and this means that they affect stars at relativelylarge radius. With the exception of the “tumbling" resonance, theremaining resonances in Tab.1 appear to play no role.In this section we identified resonances by examining the Fourierdecomposition of the perturbing potential. In first-order orbital per-turbation theory, a given Fourier term gives rise to perturbationscontaining a resonant denominator 𝜔 − 𝑛 Ω − 𝑛 Ω , where 𝜔 is thefrequency of the explicit time-dependence of the perturbing potential.If resonance is strictly defined as the vanishing of this denominator MNRAS000 , 1–17 (2020) Breen, Rozier, Heggie & Varri. (Jalali & Hunter 2005, p.810) then resonances do not occur if theperturbation grows, but “only near-resonances at which the denomi-nators [....] are small" (if the growth-rate is small). But these authorsalso (their p.819) point to a situation, closely analogous to ours, inwhich the Fourier component corresponding to the ILR is importanteven though no particle is in strict or near-resonance. It is a situa-tion in which “the denominator [...] is never large for any orbit", anapt description for the fact that the line labelled “ILR" in Fig.16 isapproximately parallel to the distribution of points. Such a Fourierterm, it seems, may have a comparable importance to a resonant termin which a small fraction of particles are in “near-resonance" but formost particles the denominator is large. Though it can be argued thatthere is no resonance at all if the denominator is not small com-pared with typical orbital frequencies, we prefer to regard the ILRas a weak resonance affecting most, if not all particles. In general, inwhat follows we take a liberal view of what is meant by “resonance".
Not all theoretical discussions of bar instabilities in disks seem read-ily adaptible to spherical systems, and here we follow specificallythe arguments proposed by Palmer, Papaloizou & Allen (1989) andAllen, Palmer & Papaloizou (1992). These two papers were con-cerned largely with two instabilities. In the earlier paper they de-veloped a theory of what they termed the circular orbit instability ,generated by the inner Lindblad resonance, and in the later paperthey turned to a so-called tumbling instability , generated by what wehave (for that reason) termed the tumbling resonance (Table 1). Bothpapers proceed from the study of the linearised, orbit-averaged col-lisionless Boltzmann - Poisson system. More precisely, they demon-strate the existence of neutral (or unstable) solutions to the resultingeigensystem by only keeping terms that correspond to a particularresonance (or a particular set of resonances) between the instabilityand the orbital frequencies in the system. We propose to use a similarmethod, and try to test whether these mechanisms can be responsiblefor instabilities in the systems studied in the present paper.
In order to investigate in detail the role of resonances in the develop-ment of instabilities in tangentially anisotropic, rotating systems, weextended the version of the matrix method formulated in Rozier, et al.(2019) to the case of rotating spheres. The matrix method identifiesunstable modes as eigenvectors (with an eigenvalue equal to 1) of theresponse matrix (cid:98) 𝑀 𝑝𝑞 ( 𝜔 ) = ( 𝜋 ) ∑︁ n ∈ Z ∫ d J n · 𝜕𝐹 / 𝜕 J 𝜔 − n · 𝛀 ( J ) (cid:2) 𝜓 ( 𝑝 ) n ( J ) (cid:3) ∗ 𝜓 ( 𝑞 ) n ( J ) , (4)where the integral runs over the action space ( 𝐽 𝑟 , 𝐿, 𝐿 𝑧 ) consisting ofthe radial action, the total angular momentum, and its 𝑧 -component,respectively; 𝐹 is the phase-space distribution function of the sys-tem; 𝛀 = ( Ω 𝑟 , Ω 𝜙 , ) are the frequencies of a given orbit; 𝜓 ( 𝑝 ) n isa Fourier-transformed element of a bi-orthogonal potential-densitybasis; n is a three-vector of integers; and 𝜔 = 𝜔 + i 𝜂 is the (complex)frequency of the instability (oscillation frequency and growth rate).Examining the integrand in Eq. (4), we can see that the denominatoris resonant, in the sense that it can take very small values in someparts of action space. Those orbits where 𝜔 − n · 𝛀 ( J ) is close to0 are expected to be decisive in the development of the instability,since they contribute the most in the response matrix. We then callthe 3-dimensional vector n of integers a resonance vector : at a given n , the corresponding integral will mainly involve the orbits that areat the n th resonance with the instability.As a consequence, each integral term in the sum of Eq. (4) can beregarded as the contribution to the response matrix from the orbitsthat are at the n th resonance with the instability. So, in order to inves-tigate the role of each resonance in the development of instabilities,one can separately study the influence of each corresponding term inthe matrix. We decided to follow this principle in order to study theinterplay between resonances and instabilities in tangentially biased,rotating systems. More specifically, instead of just comparing theamplitude of the resonance terms in Eq. (4), we opted for the fullmode search developed in Rozier, et al. (2019), applied to responsematrices in which only a subset of the resonance terms were kept:if I is a sub-ensemble of Z , we searched for modes assuming theresponse matrix is restricted to (cid:98) 𝑀 𝑝𝑞 ( 𝜔 ) = ( 𝜋 ) ∑︁ n ∈I ∫ d J n · 𝜕𝐹 / 𝜕 J 𝜔 − n · 𝛀 ( J ) (cid:2) 𝜓 ( 𝑝 ) n ( J ) (cid:3) ∗ 𝜓 ( 𝑞 ) n ( J ) . (5)Such an approach has already been used in a study of disk instabilitiesby Polyachenko (2005).An important constraint on the ensemble of resonances to considercomes from the fact that (cid:98) 𝑀 𝑝𝑞 ( 𝜔 ) ∝ 𝛿 𝑛 𝑚 𝑝 𝛿 𝑛 𝑚 𝑞 . Since we are focusingon 2-fold symmetric instabilities , we only compute matrix termswith 𝑚 𝑝 = 𝑚 𝑞 =
2, hence the sum in Eq. (5) is restricted to 𝑛 = I as asub-ensemble of Z , and define (cid:101) n = ( 𝑛 , 𝑛 ) .Using this notation, we can label some crucial resonances bytheir resonance vector (cid:101) n : (-1,2) corresponds to the inner Lind-blad resonance, (0,2) is the corotation resonance, (1,2) is the outerLindblad resonance, and so on (Table 1). Each resonance shouldbe understood through the relative value of the vector elements,so for each resonance vector (cid:101) n , its opposite − (cid:101) n corresponds tothe same kind of resonance . Hence, whenever a resonance vec-tor will be included in the ensemble I , it will be assumed thatits opposite is also included. As identified in Rozier, et al. (2019),there is little influence from | 𝑛 | > ( 𝑛 , 𝑛 ) with | 𝑛 | ≤
2. As a consequence, thecase of the complete matrix has eight different resonance vectors: I = {( , ) , ( , ) , (− , ) , ( , ) , ( , ) , (− , ) , ( , ) , ( , )} , thefirst five of which are listed in Table 1.Finally, we draw attention to the fact that a major differencebetween rotating and non-rotating systems emerges in the case (cid:101) n = ( , ) . The corresponding term in the response matrix (Eq.(4))is directly related to the 𝐿 𝑧 -gradient of the distribution function, andis independent of the orbital frequencies. We will refer to it as thetumbling term: while it does not correspond to any resonance involv-ing orbital frequencies, it is shown in Allen, Palmer & Papaloizou(1992) that this term is associated to the tumbling of orbital planesunder the effect of a perturbation (see also Palmer 1994).In summary, our aim is to investigate the role of resonances by ex-amining the effect on the growth rate and pattern speed of removingsome of the corresponding terms from the calculation. In order to Note that in principle, other instabilities could exist (i.e. with 𝑚 =
0, givingaxisymmetric oblate/prolate oscillations), however they are beyond the scopeof this paper. Note that the role of opposite resonance vectors is not the same; they onlycorrespond to the same kind of resonances. This is illustrated by the differentinclination-dependence of resonances and their “reverse" analogues in Table1.MNRAS , 1–17 (2020) tability of Rotating Anisotropic Models Table 2.
Values of the growth rate and pattern speed found through the re-stricted response matrix method applied to the ( 𝑞, 𝛼 ) = (− , . ) and ( 𝑞, 𝛼 ) = (− , ) tangentially-biased systems: comparison between the ma-trix of reference (complete) and the series of matrices obtained by removinga single resonance. I \ I (removed) 𝜂 ( . ) Ω 𝑏 ( . ) 𝜂 ( ) Ω 𝑏 ( )∅ (reference) 0.012 0.36 0.24 0.66 (− , ) ( , ) ( , ) ( , ) ( , ) (− , ) ( , ) ( , ) strengthen the case for this method, we present in Appendix C the re-sult of its application to the well-studied radial orbit instability regimefound in Rozier, et al. (2019). We now apply it to the tangentially-biased rotating systems which are the focus of the present paper, andconcentrate on the models with 𝑞 = −
16 and 𝛼 = . , In order to identify which resonances are the most critical to theinstabilities, we compared the results of the full matrix with resultswhen a single resonance is removed. The importance of the removedresonance in the instability is then estimated through the change ingrowth rate and oscillation frequency of the instabilities found withand without the resonant term.Table 2 shows the values of the growth rates and pattern speedsthat were obtained by the restricted matrix method in the two cases 𝛼 = . 𝛼 =
1, alternatively removing each of the resonantterms that are present in the matrix of reference. Notice that thevalues of the growth rate and pattern speed of reference do notexactly correspond with those of Rozier, et al. (2019): for the sakeof numerical efficiency, we opted for fewer radial basis functions,so the mode reconstruction is not as accurate. However, the presentreference case identifies the same instability, and we ensure the selfconsistency of the present method by performing matrix calculationswith constant parameters, in particular the nature and number ofradial basis functions. Similarly, the discrepancy in the measuredgrowth rate of the 𝛼 = . 𝛼 = . (− , ) , ( , ) or ( , ) does not impactthe frequencies by more than 20% . Some resonances have a moresignificant impact on the instability growth rate: removing a termamong {( , ) , ( , ) , ( , ) , ( , )} depletes the growth rate by atleast 30%, and up to 70%. The pattern speed, however, is barelyimpacted by the removal of any of these terms. It appears that thecontribution to the instability from the resonances in these groupscan be interpreted as follows: none of them has a fundamental role inthe formation of the instability, however they all contribute more or It is not hard to show that these terms are smaller than the others by a factorat least of order the epicyclic amplitude, which is small in these tangentiallyanisotropic models. less to it by increasing its growth rate. This analysis does not apply tothe ILR term (− , ) : when it is removed, the remaining instabilityhas a growth rate about 50 times smaller than the reference one,while its pattern speed is depleted by about a factor 3. This indicatesthat the instability identified in the absence of the ILR term is of adifferent nature from the one identified in all the other cases.The maximally rotating, 𝛼 = {( , ) , (− , ) , ( , ) , ( , )} .The group of resonances that influence the growth rate but not thepattern speed is now reduced to one element: the removal of the OLRterm results in the depletion of the growth rate by a factor 3 and ofthe pattern speed by about 10%. A new category emerges: the tum-bling and corotation terms are now associated to decreases of a factorrespectively about 5 and 10 in the growth rate, and of respectivelyabout 40% and 25% in the pattern speed. Finally, the ILR term is stillthe one bearing the most importance in the instability, as its removalleads to a depletion of a factor about 40 in the growth rate, and about2 in the pattern speed.The most important process that leads to instability seems, inboth 𝛼 = . 𝛼 = (− , ) , ( , ) , ( , ) andespecially ( , ) terms is lower at higher rotation, while corotation,the OLR and the tumbling terms gain importance with rotation. Thetumbling term seems to gain a particular role in setting the highpattern speed of the instability at high 𝛼 .Before we pass on from the matrix method, it is worth mentioningthat we have focused here on the most rapidly growing mode. For abroad ranges of values of large 𝛼 and negative 𝑞 (in models where theLynden-Bell demon is applied everywhere), it coexists with a secondmode which has smaller growth rate and pattern speed (by factors oforder five and two, respectively). 𝑁 -body model In this section we return to 𝑁 -body modelling to examine evidencefor resonant behaviour. In particular we discuss direct 𝑁 -body sim-ulations of basic models with 𝛾 = 𝑞 = −
16) and maximal rotation ( 𝛼 = 𝑡 =
30, but its growth is clearly nonlinear by 𝑡 =
15 (when a plot ofthe amplitude against time has an inflexion from concave upwards toconcave downwards). During this shorter interval the pattern speed Ω 𝑏 decreases from 0.73 to 0.60.The initial conditions and relevant resonances are shown in Fig.16,and now Fig.17 shows the evolution of the distribution of the az-imuthal frequency Ω during the period of roughly linear behaviour.In most cases it is difficult to distinguish features from samplingfluctuations, and we shall improve the statistics below, but in viewof the large numbers of particles per bin, it is hard to dispute thatthere is a significant increase around Ω = .
65 and an equally sig-nificant decrease around 1 .
05 (the mid-point of the histogram bar isquoted). It is tempting to link the first of these with corotation. Againtentatively, the second could be associated with the “epicyclic" reso-nance: in Fig.16 the corresponding line crosses the upper margin ofthe distribution of points, where they are concentrated, at a frequencyjust below Ω =
1. Other features at higher frequency, if regarded assignificant, correspond to none of the resonances listed in Tab.1, butthere is another at a lower frequency, about Ω = .
4, which could
MNRAS000
MNRAS000 , 1–17 (2020) Breen, Rozier, Heggie & Varri.
Figure 17.
Evolution of the distribution of Ω during the early growth of abar in an 𝑁 -body model with 𝑁 = 𝑞 = −
16 and 𝛼 =
1. Theblack bar shows the variation of Ω 𝑏 (which decreases slightly) during thesame interval of time. The green bars show the ranges of the epicyclic (right)and Outer Lindblad (left) resonances. For those resonances the breadth offrequencies is mainly due to the width of the distribution of frequencies inthe initial conditions (Fig.16). very well correspond to the outer Lindblad resonance, as can be seenfrom Fig.16.It is worth examining the distribution of those stars which con-tribute most to the evolution of 𝑓 ( Ω ) . Fig.18 shows two plots cor-responding, respectively, to low and high inclination, and also theplot for all inclinations (black). The red curve, for instance, showsthe change from 𝑡 = 𝑡 =
15 in the number of particles in abin (centred at the indicated values of Ω ), counting only those ofinclination 𝑖 < ◦ . Results from ten independent simulations havebeen combined, and the bin size is half that of Fig.17. The sum,plotted in black, adds considerable significant detail to what canbe inferred from Fig.17. In the features at Ω (cid:39) . Ω (cid:39) .
05. Bearing in mind theinclination-dependence noted in Tab.1, these results lend weight tothe interpretations as behaviour at the outer Lindblad, corotation andepicyclic resonances, respectively.Despite what has just been said about the epicyclic resonance, ac-tually the resonant frequency is somewhat too small (see the positionof the rightmost horizontal bar in the figures), and in this respecta better interpretation is offered by consideration of the ILR, theimportance of which has already been established in Sec.4.2. It isdefined by Ω 𝑏 = Ω ILR ≡ Ω − Ω /
2, and it is found that, in the limit 𝑞 → −∞ , the expression Ω ILR reaches its maximum value at a radiuswhere Ω = .
09. Its importance is emphasised again in Sec.5.1. Athird possible interpretation is that the feature near Ω = Ω , which is a decreasing functionof radius. Indeed the bins at Ω = . , .
65 and 1.05 correspondto approximate radial ranges 2 < 𝑟 < . , . < 𝑟 < . . < 𝑟 < . Ω can also be attempted with the distribution of Ω ,related to the epicyclic resonance, and Ω + Ω , involving the OLR.But again because of the narrowness of the distribution of the twobasic frequencies, one sees evidence of rapid evolution at essentiallythe same corresponding frequencies, i.e. those which may be associ- Figure 18.
Change in the distribution of Ω between 𝑡 = 𝑡 =
15, givenseparately for particles of high ( 𝑖 > ◦ ; green) and low inclination (red); thethick black line gives the sum, i.e. all inclinations. The error bars on the lattergive the average and 1- 𝜎 standard error of the mean over 10 simulations.For the first two graphs the error bars would be larger by about 40%. Thehorizontal bars at the foot have the same meaning as in Fig.17, but depend onthe bar pattern speed, which was determined for only one of the simulations.The binwidth is half that of Fig.17. ated with all three resonances, plus some evidence for unattributableevolution at still larger frequencies. This paper has presented a diversity of empirical information on barformation in a tangentially anisotropic, rotating spherical system.To tie everything together, it is helpful to set our findings within atheoretical framework, and the one which we now present, thoughnot original, has proved most fruitful in the exploratory aspects ofour research.Palmer, Papaloizou & Allen (1989) and Allen, Palmer & Pa-paloizou (1992) exhibit two particular mechanisms that can giverise to linear instability in, respectively, tangentially anisotropic androtating spherical stellar systems: the circular orbit instability andthe tumbling instability. Their approach is via a careful analysis ofthe linearised orbit-averaged Boltzmann-Poisson system, in whichthe first step is to treat one resonance in isolation, and then all othersare treated as perturbations. For both mechanisms, it is shown inthe first step that the system would exhibit a particular set of neutralmodes (i.e. modes with vanishing growth rate) in the absence of theperturbing resonances , but that such a mode may be destabilised bythe perturbing resonant interactions. While the tumbling process islinked to the precession of orbital planes under the influence of arotating bar, the circular orbit instability emerges from the influenceof a bar on the apsidal precession of the orbits, entailing changes in Ω ILR (defined in Eq.(6)).In the case of the circular orbit process, the characteristics of theinstability are the following: • The ILR is at the source of the neutral mode, hence it is thefundamental cause of the instability. • The mode rotates at a pattern speed larger than (but close to)
MNRAS , 1–17 (2020) tability of Rotating Anisotropic Models the maximum value of the ILR frequency, defined by Ω ILR = Ω − 𝜅. (6)Fig. 15 shows this frequency in the present clusters as a function ofradius. • Destabilisation is ensured by other resonances, at the oscillationfrequency of the mode.A similar list of characteristics can be given for the tumbling insta-bility: • The tumbling process (corresponding to the ( , ) resonancevector) is at the source of the neutral mode. • The pattern speed is nearly proportional to the rotation pa-rameter 𝛼 . This comes from the fact that the tumbling term in theeigenequation describing the neutral mode (see eq. (8.31) in Palmer1994) is proportional to 𝛼 / 𝜔 . • Destabilisation is ensured by other resonances, at the frequencyof the mode.Before passing on to a discussion of empirical evidence on themechanisms of bar formation, we mention here one factor whichsupports the framework of the above theory: the presence of a secondmode (see the last paragraph of Sec.4.2.2). The main problem of thecited authors’ approach, in which only one resonance is included inthe first instance, exhibits not one but a sequence of neutral modes,and there is no reason to suppose that only one of them can bedestabilised by perturbing resonances.
The series of experiments using the matrix method that we reported inSection 4.2 helps in associating the instabilities we exhibited through 𝑁 -body experiments (mainly in Sec.3) with the processes of tumblingand circular orbit instability. We begin mainly with the matrix results,and then pass on to an interpretation of the 𝑁 -body data.The ( 𝑞, 𝛼 ) = (− , . ) case is particularly eloquent : it is obvious(from Table 2) that the instability mostly depends on the ILR; thepattern speed of the mode lies just above the maximum of Ω ILR (seeFig.15); and the role of other resonances seems to be to accumulatethe growth rate. The ( 𝑞, 𝛼 ) = (− , ) case seems to follow a similartrend, yet less tightly matching the characteristics of the circularorbit instability: the ILR term still bears most of the instability;the pattern speed is twice as high as the maximum of Ω ILR ; andyet some other terms play an important role, in particular the oneassociated with the tumbling process. In contrast to the conclusionsof Allen, Palmer & Papaloizou (1992), it seems that the instabilitiesin the present tangentially anisotropic, rotating systems cannot beattributed to the tumbling instability alone. Perhaps the high valuesof the pattern speed and the growth rate may diminish the accuracyof their perturbative approach. The results of Table 2 rather point to amixing of the circular orbit and the tumbling processes, the ILR stillbeing the most important resonance in the creation of the instability.It is at first sight paradoxical to assert that the ILR is central tothe observed instabilities, when Fig.16 shows that no particles in themodels lie close to resonance. However, the restricted matrix method(Table 2) shows that, even if the pattern speed is far away from anyvalue Ω ILR can take in the system (more than twice its maximumvalue when 𝛼 = 𝑁 -body models that the ILR is playing a vital role: the frequencymeasurements presented in Figs 17 and 18 show strong activity inparticles for which Ω = .
05, which corresponds closely to the value of Ω at which Ω ILR takes its maximum value ( Ω ILR = . 𝛼 decreases; this flooris near the maximum of Ω ILR , as expected in the scenario outlinedin Sec.5.1.1.The foregoing remark, that a resonance may have a vital role evenin the absence of resonant particles, echoes a point made near the endof Sec.4.1, and should be borne in mind in assessment of the role ofthe tumbling instability. (This is formally the condition Ω 𝑏 =
0, andis never satisfied in our models). Consider in particular the region Ω > . 𝑅 < . . < 𝑅 < .
75 increases the growth rate of the instability in the innersystem up to its maximum value. By contrast, the results of these twosections also show that what happens at larger radii (
𝑅 > .
5, say)hardly matters: this region may rotate, counter-rotate, or not rotateat all. The importance of the most central orbits is also supported bythe shapes of the modes, as shown by Fig. 13 and in Rozier, et al.(2019): the modes seem to be confined to the inner part of the cluster,as by about 𝑅 = .
45 in the present units, the density in the bar hasdropped to 10% of its maximum value. All the evidence supports theview that the tumbling process, while not dominant in the generationof instability, is an important component of the mechanism whichcreates it.It should be noted that the foregoing discussion refers to the cylin-drical radius 𝑅 . For high-inclination orbits, such as those most af-fected by the tumbling and epicyclic resonances, the spherical radiusis significantly larger. Still, Fig. 18 shows that there is also someactivity in the low inclination orbits (in the inner region on which wehave been focusing).Low inclination orbits have a particularly important role in theregions of the OLR and corotation (around Ω = . 𝛼 is the theoretical result (mentioned in Sec.5.1.1) thatthe pattern speed of the relevant neutral mode is proportional to 𝛼 .This suggests an explanation of the rising pattern speed for 𝛼 > ∼ . 𝛼 , as then the dampingeffect of the ILR would come into play (Palmer, Papaloizou & Allen1989, p.1292). It is interesting to note also, from Table 2 in the case 𝛼 =
1, that removal of the ( , ) term reduces Ω 𝑏 to a value quitesimilar to the reference result for the case 𝛼 = . Ω ILR .To sum up, the present results point to the following scenariofor the formation of instabilities in tangentially anisotropic, rotatingsystems: in the inner region of the cluster, the tumbling of orbitalplanes (especially for high-inclination orbits) and the drift in apsidalprecession far from ILR cooperate to the formation of an unstable bar.The growth rate of this pattern is enhanced by resonance with outer
MNRAS , 1–17 (2020) Breen, Rozier, Heggie & Varri. orbits at the epicyclic resonance and (most noticeably) the corotationresonance and the OLR, in particular with the low-inclination orbits.Thus there is a kind of cooperation between all the resonances, whichdiminishes high-inclination orbits in the inner regions, and enhancesthose of low inclination at larger radii. Both effects create the flattenedstructure of the bar.
In their classic study Ostriker & Peebles (1973) stated that disksbecome “approximately stable" when 𝑇 /| 𝑊 | falls to a value of 0 . ± .
02. The wording shows just how difficult it is to establish stabilitythrough numerical methods, and the same problem arises in ourwork. Furthermore, the emphasis in their work was on disks (withor without halos), whereas ours is mainly (but not exclusively) onspherical models which rotate globally. Nevertheless we take theOstriker & Peebles criterion as our starting point in this subsection.
Figs. 8, 10 and 11 summarise the results of our main numericalsurveys, placing them in the space of growth rate, 𝑇 / 𝑊 and particularcircular velocities in the cluster. Even if there seems to be a trendtowards more rapid instability with increasing 𝑇 / 𝑊 , an importantfraction of our parameter space does not satisfy the Ostriker-Peeblescriterion. The present results simply add to the number of publisheddeviations from this law (e.g., see Zang & Hohl 1978; Berman &Mark 1979; Aoki, Noguchi & Iye 1979; Efstathiou et al. 1982; Evans& Read 1998; Athanassoula 2008), and yet are among the few thatexhibit such deviations in purely spherical clusters.As a consequence, it appears that the linear stability of sphericalsystems cannot be uniquely determined by the ratio between globalkinetic energy in rotating motion and global gravitational potentialenergy. Nor would the global angular momentum be more useful: oursimulations of clusters with counter-rotating components (high-shearmodels, Sec.3.2) imply the presence of fast growing instabilities inclusters with arbitrarily low global angular momentum (when theangular momentum of the inner prograde and outer retrograde partsof the cluster compensate).An important reservation should be made at this point, concerningthe presence of instabilities with low growth rates. On one hand, interms of measurements by analytical means (i.e. the response matrixin the present case), the behaviour of the analytical indicators atlow growth rate gets very intricate (see, e.g., Pichon & Cannon1997; Merritt 1999; Rozier, et al. 2019), which makes it complicatedto measure growth rates accurately. On the other, the measurementof slowly growing instabilities in numerical simulations also has anumber of technical difficulties. In particular, measuring the growthof these very weak instabilities above the background noise requiresa sufficient integration time, and in 𝑁 -body simulations as well asin actual physical systems, the time scale required for the growthof such instabilities can overlap with the secular time scale of 2-body relaxation in the cluster. As a consequence, one might miss thegrowth of an instability, because the system has already secularlyevolved to a stable configuration.Yet, even if we only focus on reasonably large growth rates (say 𝜂 > . 𝑇 / 𝑊 and the meanazimuthal velocity at 𝑅 = . 𝑇 / 𝑊 , even if the separation between stable and unstablesystems (or, rather, slowly and rapidly growing instabilities) in thatspace is not yet perfect.Let us remark in closing that the 2D space of global rotation(through 𝑇 / 𝑊 ) and azimuthal velocity at a given radius can bemapped into the space of global rotation and shear: a cluster withhigh | ¯ 𝑣 𝜙 ( . )| displays low shear if 𝑇 / 𝑊 is also high, and high shearif 𝑇 / 𝑊 is low; a cluster with low | ¯ 𝑣 𝜙 ( . )| displays high shear if 𝑇 / 𝑊 is high, and low shear if 𝑇 / 𝑊 is low. We take the view that the question of the stability of a non-rotatingequilibrium with a high degree of tangential anisotropy is still un-settled. On the one hand, Polyachenko (1987) uses the WKB ap-proximation to show that such spheres are linearly stable. However,this result only probes the case of quasi-local disturbances, while wemeasure large scale patterns (as illustrated by Fig.13), which maybe beyond the scope of the WKB approximation. On the other hand,the possibility of large scale instabilities in such spheres is supportedby Palmer, Papaloizou & Allen (1989), based on complementaryanalytical and numerical results. A consequence of this finding, ifconfirmed, would be that sufficiently tangentially anisotropic clustersseem to display unstable behaviour, whatever their degree of rota-tion or shear . But we caution that this intriguing regime deserves acareful analysis which goes beyond the scope of the current paper,and we hope to present a contribution dedicated to this topic in thefuture.
This paper focuses on the collisionless stability of rotating, spheri-cal stellar systems, using a mix of 𝑁 -body simulations, the matrixmethod for the evolution of perturbations, and some physical ideas.Our starting point is a set of anisotropic Plummer models (Dejonghe1987), whose anisotropy is parametrised by a parameter 𝑞 ; our focusis on models with 𝑞 <
0, whose velocity distribution is tangentiallybiased. To these models we apply different versions of the Lynden-Bell demon (Lynden-Bell 1960), in which the rotation of a fraction 𝛼 of stars is reversed in a region of phase space defined in terms of inte-grals of motion, and parametrised in various ways by a parameter 𝛾 .Instances of these models were generated with the publicly availablesoftware PlummerPlus , and their evolution studied using the publicversion of gyrfalcON . Supplementary studies used NBODY6 (Nita-dori & Aarseth 2012). The matrix method, described for example inRozier, et al. (2019), was also used to search for unstable modes.Our conclusions are as follows:(i) For values of 𝑞 < < 𝛼 ≤ 𝑞 = −
6, and afraction 𝛼 of retrograde stars have the rotation reversed, the growthrate of the mode is approximately proportional to 𝛼 (Fig.5), atleast when 𝛼 is large enough that the growth rate can be measuredreliably. In extremely tangentially anisotropic models, the pattern https://teuben.github.io/nemo/man_html/gyrfalcON.1.htmlMNRAS , 1–17 (2020) tability of Rotating Anisotropic Models speed diminishes as 𝛼 decreases, but remains above a certain limit,and is detectable even in non-rotating systems (Fig.14).(iii) If only stars with energy below a certain fraction 𝛾 of theminimum energy have their rotation reversed, the growth rate reachesroughly its maximum value if 𝛾 > ∼ .
5, i.e. if the Lynden-Bell demonis applied to at least the most bound half of the stars (Fig.6).(iv) In models in which stars below a given energy have theLynden-Bell demon applied in a prograde sense and for the remainderit is applied in a retrograde sense, the minimum growth rate occurswhen the proportions of retrograde and prograde stars in the mostbound 40% of the system are nearly the same (i.e. the case 𝛾 = . 𝑅 < ∼ . , | 𝑧 | < . rapid instability, we showthat, in addition to the global rotational kinetic energy, it is helpfulto include a measure of the rotation speed at a radius inside the virialradius of the system (Fig.11). ACKNOWLEDGEMENTS
We thank an anonymous referee for a number of helpful com-ments, including pertinent references. SR thanks Christophe Pichonand Jean-Baptiste Fouvry for many interesting discussions, and allauthors are grateful to them for commenting on an early draft.This work has made use of the Horizon Cluster hosted by Institutd’Astrophysique de Paris. All authors thank Stephane Rouberol forrunning smoothly this cluster. We are grateful to S. Law and his ITstaff in the Mathematics Department at Edinburgh University, forvery helpful management of the GPU hardware. Many simulationswere also carried out with the ECDF cluster at the University. Allauthors thank the Leverhulme Trust for financial support from grantRPG-2015-408. ALV acknowledges support from a UKRI FutureLeaders Fellowship (MR/S018859/1) and is grateful to the Depart-ment of Astronomy at the University of Tokyo, where part of thiswork was conducted, thanks to a JSPS International Fellowship withGrant-in-Aid (KAKENHI-18F18787).
Data Availability
Data underlying this article will be shared on reasonable requestto the authors via [email protected]. Data related to the initialconditions may be reproduced via the publicly available software
PlummerPlus . REFERENCES
Alimi J.-M., Perez J., Serna A., 1999, MNRAS, 305, 859Allen A. J., Palmer P. L., Papaloizou J. C. B., 1992, MNRAS, 256, 695Anderson J., King I. R., 2003, AJ, 126, 772Antonov V. A., 1987, IAUS, 127, 531Aoki S., Noguchi M., Iye M., 1979, PASJ, 31, 737Athanassoula, E., 2008, MNRAS, 390, L69Barnes J., Goodman J., Hut P., 1986, ApJ, 300, 112Bellini A., et al., 2017, ApJ, 844, 167Berman, R. H., Mark, J. W.-K., 1979, A&A, 77, 31Bertin G., Pegoraro F., Rubini F., Vesperini E., 1994, ApJ, 434, 94Bertin, G., Lin, C. C., 1996, Spiral Structure in Galaxies: a Density WaveTheory, (Cambridge, MA: MIT Press)Bianchini P., et al., 2018, MNRAS, 481, 2125Binney J., Tremaine S., 2008, Galactic Dynamics, 2e (Princeton:PrincetonUniversity Press)Breen P. G., Varri A. L., Heggie D. C., 2017, MNRAS, 471, 2778Centrella, J. M., New, K. C., Lowe, L. L., et al., 2001, ApJ, 550, L193Chandrasekhar, S., 1969, Ellipsoidal Figures of Equilibrium (New Haven,CT: Yale Univ. Press)Christodoulou, D. M., Shlosman, I., Tohline, J. E., 1995, ApJ, 443, 551Dehnen W., 2002, JCoPh, 179, 27Dejonghe H., 1987, MNRAS, 224, 13Dejonghe H., Merritt D., 1988, ApJ, 328, 93Doremus J.-P., Feix M. R., Baumann G., 1971, PhRvL, 26, 725Efstathiou, G., Lake, G., Negroponte, J., 1982, MNRAS, 199, 1069Evans N. W., Read J. C. A., 1998, MNRAS, 300, 106Ferraro F. R., et al., 2018, ApJ, 860, 50Fujiwara T., 1983, PASJ, 35, 547Gaia Collaboration, et al., 2018, A&A, 616, A12Gillon D., 1980, CRASB, 290, 545Hachisu, I., 1986, ApJS, 61, 479Heggie D. C., Breen P. G., Varri A. L., 2020, MNRAS, 492, 6019Hamilton C., Fouvry J.-B., Binney J., Pichon C., 2018, MNRAS, 481, 2041Hénon M., 1973, A&A, 24, 229Hernquist L., Ostriker J. P., 1992, ApJ, 386, 375Hohl, F., 1971, ApJ, 168, 343Hohl, F., 1976, AJ, 81, 30Jalali M. A., Hunter C., 2005, ApJ, 630, 804James, R. A., 1964, ApJ, 140, 552Jindal A., Webb J. J., Bovy J., 2019, MNRAS, 487, 3693Kalnajs, A. J., 1972, ApJ, 175, 63Kalnajs A. J., 1977, ApJ, 212, 637Kamann S., et al., 2018, MNRAS, 473, 5591Karino, S., Eriguchi, Y., 2003, ApJ, 592, 1119Kuijken, K., Dubinski, J., 1994, MNRAS, 269, 13Lai, D., Rasio, F. A., Shapiro, S. L., 1993, ApJS, 88, 205Lin, C. C., 1955, The Theory of Hydrodynamic Stability, (Cambridge, UK:Cambridge Univ. Press)Lynden-Bell D., 1960, MNRAS, 120, 204Lynden-Bell D., 1962, MNRAS, 123, 447Lynden-Bell D., 1969, MNRAS, 144, 189Lynden-Bell D., Sanitt N., 1969, MNRAS, 143, 167Lynden-Bell D., Kalnajs A. J., 1972, MNRAS, 157, 1Mark J. W.-K., 1977, ApJ, 212, 645Merritt D., 1999, PASP, 111, 129Merritt D., Aguilar L. A., 1985, MNRAS, 217, 787Meza A., 2002, A&A, 395, 25 MNRAS000
Alimi J.-M., Perez J., Serna A., 1999, MNRAS, 305, 859Allen A. J., Palmer P. L., Papaloizou J. C. B., 1992, MNRAS, 256, 695Anderson J., King I. R., 2003, AJ, 126, 772Antonov V. A., 1987, IAUS, 127, 531Aoki S., Noguchi M., Iye M., 1979, PASJ, 31, 737Athanassoula, E., 2008, MNRAS, 390, L69Barnes J., Goodman J., Hut P., 1986, ApJ, 300, 112Bellini A., et al., 2017, ApJ, 844, 167Berman, R. H., Mark, J. W.-K., 1979, A&A, 77, 31Bertin G., Pegoraro F., Rubini F., Vesperini E., 1994, ApJ, 434, 94Bertin, G., Lin, C. C., 1996, Spiral Structure in Galaxies: a Density WaveTheory, (Cambridge, MA: MIT Press)Bianchini P., et al., 2018, MNRAS, 481, 2125Binney J., Tremaine S., 2008, Galactic Dynamics, 2e (Princeton:PrincetonUniversity Press)Breen P. G., Varri A. L., Heggie D. C., 2017, MNRAS, 471, 2778Centrella, J. M., New, K. C., Lowe, L. L., et al., 2001, ApJ, 550, L193Chandrasekhar, S., 1969, Ellipsoidal Figures of Equilibrium (New Haven,CT: Yale Univ. Press)Christodoulou, D. M., Shlosman, I., Tohline, J. E., 1995, ApJ, 443, 551Dehnen W., 2002, JCoPh, 179, 27Dejonghe H., 1987, MNRAS, 224, 13Dejonghe H., Merritt D., 1988, ApJ, 328, 93Doremus J.-P., Feix M. R., Baumann G., 1971, PhRvL, 26, 725Efstathiou, G., Lake, G., Negroponte, J., 1982, MNRAS, 199, 1069Evans N. W., Read J. C. A., 1998, MNRAS, 300, 106Ferraro F. R., et al., 2018, ApJ, 860, 50Fujiwara T., 1983, PASJ, 35, 547Gaia Collaboration, et al., 2018, A&A, 616, A12Gillon D., 1980, CRASB, 290, 545Hachisu, I., 1986, ApJS, 61, 479Heggie D. C., Breen P. G., Varri A. L., 2020, MNRAS, 492, 6019Hamilton C., Fouvry J.-B., Binney J., Pichon C., 2018, MNRAS, 481, 2041Hénon M., 1973, A&A, 24, 229Hernquist L., Ostriker J. P., 1992, ApJ, 386, 375Hohl, F., 1971, ApJ, 168, 343Hohl, F., 1976, AJ, 81, 30Jalali M. A., Hunter C., 2005, ApJ, 630, 804James, R. A., 1964, ApJ, 140, 552Jindal A., Webb J. J., Bovy J., 2019, MNRAS, 487, 3693Kalnajs, A. J., 1972, ApJ, 175, 63Kalnajs A. J., 1977, ApJ, 212, 637Kamann S., et al., 2018, MNRAS, 473, 5591Karino, S., Eriguchi, Y., 2003, ApJ, 592, 1119Kuijken, K., Dubinski, J., 1994, MNRAS, 269, 13Lai, D., Rasio, F. A., Shapiro, S. L., 1993, ApJS, 88, 205Lin, C. C., 1955, The Theory of Hydrodynamic Stability, (Cambridge, UK:Cambridge Univ. Press)Lynden-Bell D., 1960, MNRAS, 120, 204Lynden-Bell D., 1962, MNRAS, 123, 447Lynden-Bell D., 1969, MNRAS, 144, 189Lynden-Bell D., Sanitt N., 1969, MNRAS, 143, 167Lynden-Bell D., Kalnajs A. J., 1972, MNRAS, 157, 1Mark J. W.-K., 1977, ApJ, 212, 645Merritt D., 1999, PASP, 111, 129Merritt D., Aguilar L. A., 1985, MNRAS, 217, 787Meza A., 2002, A&A, 395, 25 MNRAS000 , 1–17 (2020) Breen, Rozier, Heggie & Varri.
Milone A. P., Marino A. F., Mastrobuono-Battisti A., Lagioia E. P., 2018,MNRAS, 479, 5005Mikhaˇilovskiˇi A. B., Fridman A. M., Epel’Baum Y. G., 1971, JETP, 32, 878Morozov A. G., Polyachenko V. L., Shukhman I. G., 1974, SvA, 18, 44New, K. C., Centrella, J. M., Tohline, J. E., 2000, Phys. Rev. D, 62, 064019Nitadori K., Aarseth S. J., 2012, MNRAS, 424, 545Orr, W. McF., 1907, Proc. Roy. Irish Acad. Sect. A, 27, 9Ostriker, J. P., Mark, J. W. -K., 1968, ApJ, 151, 1075Ostriker J. P., Peebles P. J. E., 1973, ApJ, 186, 467Ou, S., Tohline, J. E., 2006, ApJ, 651, 1068Palchik M. I., Patashinskii A. Z., Pinus V. K., 1974, PMTF, 63Palmer P. L., 1994, Stability of collisionless stellar systems (Dordrecht, NL:Kluwer Acad. Publishers, ASSL)Palmer P. L., Papaloizou J., 1987, MNRAS, 224, 1043Palmer P. L., Papaloizou J., Allen A. J., 1989, MNRAS, 238, 1281Palmer P. L., Papaloizou J., Allen A. J., 1990, MNRAS, 243, 282Perez J., Alimi J.-M., Aly J.-J., Scholl H., 1996, MNRAS, 280, 700Perez J., Aly J.-J., 1996, MNRAS, 280, 689Pichon C., Cannon R. C., 1997, MNRAS, 291, 616Pickett, B. K., Durisen, R. H., Davis, G. A., 1996, ApJ, 458, 714Polyachenko V. L., 1987, IAUS, 127, 301Polyachenko V. L., Shukhman I. G., 1981, SvA, 25, 533Polyachenko E. V., 2005, MNRAS, 357, 559Polyachenko E. V., Polyachenko V. L., Shukhman I. G., 2007, MNRAS, 379,573Polyachenko E. V., Shukhman I. G., 2015, MNRAS, 451, 601Rozier S., 2020, PhD ThesisRozier S., Fouvry J.-B., Breen P. G., Varri A. L., Pichon C., Heggie D. C.,2019, MNRAS, 487, 711Saha P., 1991, MNRAS, 248, 494Saijo, M., Baumgarte, T. W., Shapiro, S. L., 2003, ApJ, 595, 352Saijo, M., Yoshida, S., 2006, MNRAS, 368, 1429Saijo M., Yoshida S., 2016, PhRvD, 94, 084032Sellwood, J. A., 1981, A&A, 99, 362Sellwood J. A., 1997, ASPC, 123, 215Sellwood, J. A., Valluri, M., 1997, MNRAS, 287, 124Shibata, M., Karino, S., Eriguchi, Y., 2002, MNRAS, 334, L27Synakh V. S., Fridman A. M., Shukhman I. G., 1971, DoSSR, 201, 827Sollima A., Baumgardt H., Hilker M., 2019, MNRAS, 485, 1460Stoeckly, R., 1965, ApJ, 142, 208Tassoul J.-L., 1978, Theory of Rotating Stars, (Princeton: Princeton Univ.Press)Tohline, J. E., Durisen, R. H., McCollough, M., 1985, ApJ, 298, 220Toomre A., 1981, in S.M. Fall and D. Lynden-Bell (Eds), The Structure andEvolution of Normal Galaxies (Cambridge, UK: Cambridge Univ. Press),111Vandervoort, P. O., 1980, ApJ, 240, 478Vandervoort, P. O., 1982, ApJ, 256, L41Vandervoort, P. O., Welty, D. E., 1982, ApJ, 263, 654Vandervoort, P. O., 1991, ApJ, 377, 49Watts, A. L., Andersson, N., Beyer, H., et al., 2003, MNRAS, 342, 1156Watts, A. L., Andersson, N., Williams, R. L., 2004, MNRAS, 350, 927Watts, A. L., Andersson, N., Jones, D. I., 2005, ApJ, 618, L37Weinberg M. D., 1994, ApJ, 421, 481Yoshida S., Saijo M., 2017, MNRAS, 466, 600Zang, T. A., Hohl, F., 1978, ApJ, 226, 521
APPENDIX A: DETERMINATION OF TIME SCALES INTHE MODELS OF SEC.3.4.2
This appendix describes details of the way in which data for Fig.14were derived from results of 𝑁 -body simulations.The analysis of 𝐶 (Eq.(2)) began with the growth rate. For large 𝛼 , | 𝐶 𝑚 | was found to saturate at some time 𝑡 𝑚𝑎𝑥 . For these models( 𝛼 ≥ .
6) the growth rate was obtained by fitting a logistic curve(which is nearly exponential at early times) in this time interval. For smaller 𝛼 , where the growth of a bar (if present) did not saturateby time 𝑡 = 𝛼 = . 𝛼 = , .
2. In all cases, error bars represent the 1- 𝜎 asymptotic confidence interval.The pattern speed was estimated by examining the time series ofarg ( 𝐶 ) . For a rotating bar, this exhibits a sawtooth wave form, theangular speed of each segment being twice the pattern speed. Forlarge 𝛼 it was evaluated by simply counting the number of periodswithin the selected time interval (see above), and the motion couldbe easily seen to be prograde. For 𝛼 = . ( 𝐶 ) . Averaging resultsover several runs (as mentioned above) gave a reliable autocorrela-tion with an estimate of the standard error. The form of the averageautocorrelation function is a clear decaying oscillation. The locationof its first maximum was determined by applying a quadratic fit tothe section on either side of the maximum, as far as the neighbouringminima. This method does not allow determination of the sense ofpattern rotation. This is therefore undetermined for 𝛼 = .
2, but ispresumed to be prograde. For 𝛼 = 𝑙 =
2, but three are filtered out by the definition of 𝐶 ). APPENDIX B: RESONANCES IN INCLINED EPICYCLICMOTIONS PERTURBED BY A ROTATING BAR
In the context of disk dynamics, Binney & Tremaine (2008, Sec.3.3.3)use the expression Φ ( 𝑅, 𝜑 ) = Φ 𝑏 ( 𝑅 ) cos [ ( 𝜑 − Ω 𝑏 𝑡 )] as a modelfor the potential of a planar bar, where 𝑅, 𝜑 are plane polar coordi-nates in the orbital plane, Ω 𝑏 is the pattern speed, and Φ 𝑏 specifiesthe radial dependence of the potential. The obvious extension tothree dimensions is Φ ( 𝑟, 𝜃, 𝜙 ) = Φ 𝑏 ( 𝑟 ) sin 𝜃 cos [ ( 𝜙 − Ω 𝑏 𝑡 )] inspherical coordinates, since the angle-dependence is that of an 𝑙 = 𝑖 to the 𝑥, 𝑦 plane ( 𝜃 = 𝜋 / Ω at time 𝑡 = . In cylindrical polar coordinates ( 𝑅, 𝜑, 𝑧 ) based on the orbital plane, with the angular coordinate 𝜑 = Φ ( 𝑅, 𝜑, 𝑧 = ) = Φ 𝑏 ( 𝑅 ) (cid:110) 𝑖 cos [ ( Ω − Ω 𝑏 𝑡 )]+( + cos 𝑖 ) cos ([ ( Ω − Ω 𝑏 𝑡 + 𝜑 )] + (B1) ( − cos 𝑖 ) cos ([ ( Ω − Ω 𝑏 𝑡 − 𝜑 )] (cid:111) . Now we make use of the epicyclic approximation, by which (Bin-ney & Tremaine 2008, eqs.(3.91),(3.93a)) 𝑅 = 𝑅 + 𝑋 cos 𝜃 (B2) 𝜑 = 𝜃 − 𝛾𝑋 sin 𝜃 , (B3)where ( 𝑅 , 𝜃 ) are the cylindrical radius and longitude of the guidingcentre, 𝑋 is the radial amplitude of the epicyclic motion, 𝜃 is the In the rest of the paper Ω means the circular angular velocity.MNRAS , 1–17 (2020) tability of Rotating Anisotropic Models Table C1.
Values of the growth rate and pattern speed found through thetruncated response matrix method applied to the 𝑞 = , 𝛼 = I 𝜂 Ω 𝑝 Reference 0.088 0.073 (− , ) , ( , ) I \ (− , ) < × − - I \ ( , ) (− , ) ( , ) < × − - phase of the epicycle, and 𝛾 is a certain function of 𝑅 . SubstitutingEqs.(B2) and (B3) into Eq.(B1) and expanding to lowest order in 𝑋 ,we can readily locate all trigonometrical terms of the form (cid:26) cossin (cid:27) ( ( Ω − Ω 𝑏 𝑡 ) + 𝑛 𝜃 + 𝑛 𝜃 ) , (B4)where 𝑛 , 𝑛 are (small) integers. This includes the pair (0,0), whicharises from the first term in Eq.(B1). Using the fact that (cid:164) 𝜃 , =Ω , (the radial and azimuthal frequencies), these terms yield theresonances listed in Table 1, including the stated 𝑖 -dependence andthe values in the last column. APPENDIX C: APPLICATION OF THE RESTRICTEDMATRIX TO THE RADIAL ORBIT INSTABILITY
As described in Binney & Tremaine (2008), spherical systems bi-ased towards radial orbits are subject to the so-called radial orbitinstability. This instability emerges through a Jeans process, by theclustering of radially elongated orbits around a bar-like overdensity.More precisely, the torque applied by a weak bar on a radially elon-gated orbit acts to modify the apsidal precession rate of the orbit (i.e. Ω − Ω / 𝑞 = , 𝛼 = | 𝑛 | , | 𝑛 | ≤
2, correspondingto I = {( , ) , ( , ) , ( , ) , (− , ) , (− , ) , ( , ) , ( , ) , ( , )} . Itappears that turning off terms corresponding either to the ILR or tothe tumbling process from the full matrix lowers the growth rateand modifies the oscillation frequency of the instability, while keep-ing only those two terms together preserves both its growth rateand (nearly) its oscillation frequency. These results show that thedestabilisation of radially anisotropic rotating systems is guaranteedthrough the interaction of the radial orbit and tumbling processesonly. In that case, the identification of these processes is made easyby the largely dominant role played by the corresponding terms inthe response matrix: the second line of Table C1 shows that all otherresonance vectors can be neglected w.r.t. (− , ) , ( , ) . Section 4.2 shows that tangentially anisotropic systems are more complex in thatrespect. This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000