# The kinematic richness of star clusters – II. Stability of spherical anisotropic models with rotation

Philip Breen, Simon Rozier, Douglas Heggie, Anna Lisa Varri

► To cite this version:

Philip Breen, Simon Rozier, Douglas Heggie, Anna Lisa Varri. The kinematic richness of star clusters – II. Stability of spherical anisotropic models with rotation. Monthly Notices of the Royal Astronomical Society, 2021, 502 (4), pp.4762-4778. 10.1093/mnras/stab365. hal-03176045

HAL Id: hal-03176045

https://hal.science/hal-03176045

Submitted on 12 Aug 2022

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

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

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

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

# The kinematic richness of star clusters – II. Stability of spherical anisotropic models with rotation

Philip G. Breen,<fnref n="1" /> Simon Rozier,<fnref n="2" /> Douglas C. Heggie and Anna Lisa Varri<fnref n="1" />,<fnref n="3" />

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

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 *N*-body simulations and also the matrix technique for calculation of modes with the perturbed collisionless Boltzmann equation. These methods are applied to spherical systems with an initial Plummer density distribution, but modified kinematically in two ways: the velocity distribution is tangentially anisotropic, using results of Dejonghe, and the system is set in rotation by reversing the velocities of a fraction of stars in various regions of phase space, à la Lynden-Bell. The aim of the *N*-body simulations is first to survey the parameter space, and, using those results, to identify regions of phase space (by radius and orbital inclination) that have the most important influence on the bar instability. The matrix method is then used to identify the resonant interactions in the system that have the greatest effect on the growth rate of a bar. Complementary series of *N*-body simulations examine these processes in relation to the evolving frequency distribution and the pattern speed. Finally, the results are synthesized with an existing theoretical framework, and used to consider the old question of constructing a stability criterion.

**Key words:** galaxies: kinematics and dynamics – galaxies: star clusters: general.

## 1 INTRODUCTION

The analysis of the stability properties of self-gravitating rotating equilibria is a classical problem in fluid and stellar dynamics (Chandrasekhar 1969), starting with the study of the ellipsoidal figures of equilibrium. The sequence of the uniformly rotating axisymmetric Maclaurin spheroids, as parametrized by the ratio of the rotational kinetic energy to the gravitational energy $T/|W|$, exhibits a bifurcation point (at $T/|W| = 0.13572$), where a new sequence of triaxial solutions branches off (Jacobi ellipsoids). Maclaurin spheroids with values of $T/|W|$ that are greater than such a threshold are *secularly* unstable with respect to global bar modes. Further along the Maclaurin sequence (when $T/|W| > 0.2738$), a *dynamical* instability with respect to a bar mode sets in.

In the theory of rotating stars (see e.g. Tassoul 1978 for a comprehensive introduction), the investigation of the properties of self-gravitating rotating fluid bodies has been generalized to the case of configurations with non-uniform density, especially polytropic fluids with solid-body (see e.g. James 1964; Stoeckly 1965; Lai, Rasio & Shapiro 1993) and with differential rotation (see e.g. Ostriker & Mark 1968; Tohline, Durisen & McCollough 1985; Hachisu 1986; Pickett, Durisen & Davis 1996; New, Centrella & Tohline 2000).

On the other hand, in the study of rotating stellar dynamical systems, stability analyses have been limited mainly to the context of discs. Hohl (1971) and Kalnajs (1972) provided the numerical and analytical evidence that uniformly rotating discs can be strongly

unstable with respect to global bar modes, if the ordered kinetic energy dominates the energy budget of the system. Ostriker & Peebles (1973) extended this numerical investigation to the case of differentially rotating discs and conjectured the condition $T/|W| < 0.14 \pm 0.02$ as an *empirical* necessary (but not sufficient) criterion for the dynamical stability of any rotating stellar system with respect to global bar-like modes.

Over the years, a number of studies, mostly based on *N*-body simulations, have been carried out to investigate the stability of stellar discs (see e.g. Hohl 1976; Sellwood 1981), so that the nature of the rotational instability with respect to bar-like modes was indeed confirmed to be dynamical. Few attempts at providing a physical interpretation of the Ostriker & Peebles (1973) criterion have been made (Vandervoort 1982); several counterexamples or alternative stability criteria have also been suggested (see e.g. Zang & Hohl 1978; Aoki, Noguchi & Iye 1979; Berman & Mark 1979; Efstathiou, Lake & Negroponte 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) and the connection with the corresponding fluid systems is only partially understood (see Christodoulou, Shlosman & Tohline 1995). Some exceptions are represented by the collisionless counterparts of the uniformly rotating polytropes (Vandervoort 1980), Riemann ellipsoids (Vandervoort & Welty 1982), and Maclaurin spheroids (Vandervoort 1991), for which analytical studies have demonstrated that the dynamical instability with respect to bar mode sets in at lower values of rotation with respect to the Maclaurin fluid sequence (e.g. at $T/|W| = 0.17114$, see Christodoulou et al. 1995).

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

<footnote n="2">²Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France</footnote>

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

<footnote>* E-mail: d.c.heggie@ed.ac.uk</footnote>
More recently, there has been a revival of interest in the study of the stability of differentially rotating fluids, kindled by the surprising discovery by Centrella et al. (2001) of an unstable $m = 1$ mode in polytropes with strong differential rotation, polytropic index $n = 3.33$, and $T/|W| \approx 0.14$. Since then, numerical studies have confirmed that $m = 1, 2$ modes can become unstable in a variety of differentially rotating fluid models, having values of $T/|W|$ as low as 0.01 (Shibata, Karino & Eriguchi 2002; Karino & Eriguchi 2003; Saijo, Baumgarte & Shapiro 2003). The study of the stability of differentially rotating spherical shells (Watts et al. 2003; Watts, Andersson & Williams 2004) suggests that the unstable modes within a low-$T/|W|$ differentially rotating configuration are characterized by corotation within the system. In particular, Watts, Andersson & Jones (2005) argued that differential rotation must be sufficiently strong in order for instability to occur. Numerical investigations of the stability of differentially rotating polytropes by Saijo & Yoshida (2006) and Ou & Tohline (2006) strengthened this picture, which has now been formalized also in an analytical setting (Saijo & Yoshida 2016; Yoshida & Saijo 2017).

In a broader stellar dynamical context, it has been known for over a century that corotation plays a central role in the stability of shear flows (see e.g. Orr 1907; Lin 1955); in the case of self-gravitating stellar discs, this intuition has offered a key guiding principle since the beginning of the development of the density wave theory of spiral structure (see e.g. Bertin & Lin 1996; Binney & Tremaine 2008, and references therein), especially in the context of various possible mechanisms underpinning the formation of a central bar in disc galaxies (see e.g. pioneering studies by Lynden-Bell & Kalnajs 1972; 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 investigations have been focused on the behaviour of spherical, non-rotating stellar systems in a purely collisionless setting. Also in this case, the analogy with simple fluid systems has offered a fruitful pathway: a numerical stability analysis of the collisionless counterparts of classic polytropes has been conducted by Hénon (1973), while an analytical criterion for the stability of the case of polytropic index $n = 3/2$ has been proposed by Antonov (1987). Such studies are also representative examples of the two main approaches broadly adopted in this context: direct numerical integration of the relevant equations of motion or of the corresponding collisionless Boltzmann equation (see e.g. Barnes, Goodman & Hut 1986; Fujiwara 1983, respectively) and variational principles in energy space (see e.g. Lynden-Bell 1969; Lynden-Bell & Sanitt 1969; Doremus, Feix & Baumann 1971; Gillon 1980). More recently, the stability analysis of isotropic spheres has also been extended beyond the search for bar modes, and the possibility of the existence of lopsided, weakly decaying modes has been revealed (see Weinberg 1994; Heggie, Breen & Varri 2020, and references therein).

Subsequent investigations of the behaviour of collisionless spherical equilibria have progressively ventured into the anisotropic regime, mostly by means of techniques based on the analysis of the relevant response matrix or a variety of N-body approaches. Much attention has been devoted to the case of radially biased systems (see e.g. Merritt & Aguilar 1985; Palmer & Papaloizou 1987; Dejonghe & Merritt 1988; Saha 1991; Bertin et al. 1994; Polyachenko, Polyachenko & Shukhman 2007; Polyachenko & Shukhman 2015 among many others), possibly under the stimulus offered by the identification of a global criterion for the so-called radial orbit instability, first proposed by Polyachenko & Shukhman (1981). On the other hand, tangentially anisotropic systems have enjoyed a lower degree of attention, with only a few studies addressing the non-

trivial investigation of their stability (see especially Mikhaǐlovskiǐ, Fridman & Epel'Baum 1971; Palchik, Patashinskii & Pinus 1974; Barnes et al. 1986; Polyachenko 1987; Palmer, Papaloizou & Allen 1989). Selected studies have also considered the case of spherically symmetric, rotating, collisionless systems, primarily with approaches based on the linear response theory (see e.g. the early works on homogeneous systems with rotation by Synakh, Fridman & Shukhman 1971; Morozov, Polyachenko & Shukhman 1974) or the direct numerical integration of the equations of motion (see e.g. Perez et al. 1996; Perez & Aly 1996; Alimi, Perez & Serna 1999; Meza 2002). A comprehensive description of several results achieved in this area of research is offered by Polyachenko & Shukhman (1981) and Palmer (1994). For the purpose of the analysis conducted here, we wish to highlight in particular two mechanisms that can give rise to linear instability in, respectively, tangentially anisotropic and rotating spherical stellar systems: the 'circular orbit instability' (Palmer et al. 1989) and the 'tumbling instability' (Allen, Palmer & Papaloizou 1992).

This study offers a further step along this classic line of theoretical development, by conducting a stability analysis of several classes of equilibrium models characterized by spherical symmetry, anisotropy in the velocity space and differential rotation, or shear. The properties of the equilibria under consideration certainly cannot capture, in full generality, the problem of investigating the stability of spheroidal, anisotropic, rotating stellar systems, none the less they offer a useful setting to bridge our fundamental understanding of axisymmetric systems with different shear prescriptions and spherical systems with various degrees and flavours of anisotropy. Particular attention is paid to the case of tangentially anisotropic systems, which, compared to the radially anisotropic ones, have been explored to a much lesser extent. The analysis of the tangentially anisotropic regime has certainly much to reveal yet, as evidenced by a new sequence of globally unstable, growing modes that has been recently identified by Rozier et al. (2019).

This study is also part of a series of articles devoted to the study of 'kinematic richness' (i.e. any deviations from the simplifying assumptions of isotropy in the velocity space and absence of internal rotation), which is empirically emerging in a variety of stellar systems, including star clusters. Indeed, the second Data Release of the European Space Agency mission Gaia has opened the era of 'precision astrometry' for the study of Galactic astronomy (Gaia Collaboration 2018). This new generation of astrometric data, complemented by proper motion measurements from decades-long campaigns with the Hubble Space Telescope (see e.g. Anderson & King 2003; Bellini et al. 2017) and state-of-the-art spectroscopic surveys (see e.g. Ferraro et al. 2018; Kamann et al. 2018), enable us to unlock, for the first time, the full 'phase space' of Galactic globular star clusters in the 6D position and velocity space. Such recent kinematic studies have also further confirmed the growing evidence that the presence of internal rotation in globular clusters is much more common than previously assumed (see e.g. Bianchini et al. 2018; Sollima, Baumgardt & Hilker 2019) and have enabled direct measurements of the degree of anisotropy of their velocity space (Milone et al. 2018; Jindal, Webb & Bovy 2019). This rapidly evolving observational landscape therefore offers the motivation to address a number of old and new questions concerning the stability properties of this class of stellar systems, with a focus on the role of internal rotation and velocity anisotropy. More specifically, this contribution generalizes the spherical equilibria adopted by Breen, Varri & Heggie (2017) for the study of the role of primordial anisotropy in the dynamical evolution of collisional stellar systems (Part I) and prepares the ground for a further evolutionary analysis
that will enrich our current fundamental understanding of the role of angular momentum in collisional stellar dynamics (Part III; Breen et al., in preparation).

The present article is organized as follows. Section 2 describes the classes of equilibria adopted as initial conditions of our *N*-body simulations and outlines the methodology employed for the study and characterization of the relevant unstable modes. The results of our stability analysis are presented in Section 3, together with a study of the spatial-temporal structure of the instability. Next, in Section 4, we investigate the role of resonances, supplementing the analysis of our *N*-body models with a comparative study conducted with the aid of the matrix method (as customized to rotating spheres by Rozier et al. 2019). Finally, in Section 5, we discuss our results within existing theoretical mechanisms for bar formation and the broader context of global stability criteria for rotating, anisotropic stellar systems.

# 2 METHOD

## 2.1 Matrix and *N*-body methods

For the study of collisionless stability in stellar systems, the two most important methods are *N*-body techniques and matrix methods. The latter deal only with linear stability, and tend to be less flexible, but provide a direct insight into the dynamics, being naturally expressible in terms of action-angle variables (or other invariants). We therefore adopt this approach in Section 4.2, where we aim to study the role of various resonances in the observed instabilities. In order to explore a wide range of initial conditions (especially in the rotation profile), however, the results we present first (Sections 3.1–3.3) adopt particle methods. The following paragraphs consider these two kinds of techniques in turn.

Generally, we avoid direct-summation *N*-body methods, because the need to complete many simulations in a reasonable time tends to restrict systems to modest *N*, when stability may be modified by collisional processes. It has even been said that they ‘should be regarded as methods of last resort’ (Sellwood 1997). Therefore, we mainly adopt an efficient (hierarchical) method which makes it feasible to carry out simulations with larger particle numbers, and employs softening to further reduce collisionality; for these reasons our method of first resort in this work is gyrfalcON (Dehnen 2002). This is used to generate the main bulk of the numerical results, which are in the earlier parts of Section 3. Elsewhere, however, in Sections 3.4 and 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 (though his emphasis was on discs) and to the study of non-rotating spherical systems by Polyachenko & Shukhman (1981). For such systems, the formulation we employ (Rozier et al. 2019; see, for example, equation 4 in this paper) would be identical to that in the appendix of Polyachenko & Shukhman (1981), except for notation, which in our case is derived from Hamilton et al. (2018). Our choice of basis functions (of biorthogonal density-potential pairs) can be found in Hernquist & Ostriker (1992) and relies on ultraspherical polynomials. In the case of rotating systems, our formulation of the matrix equation may also differ a little from those found in the literature (e.g. Palmer 1994, equation 8.18), not only because of notation, but also through the choice of distribution function, though such differences are minor.

Throughout the paper (even in sections not devoted to *N*-body simulations), units are such that *G*, the total mass and virial radius are unity (i.e. Hénon Units).

## 2.2 Initial conditions

The initial conditions we adopt all have the same density profile, that of a Plummer model, and the velocity distribution is constructed in a series of steps:

(i) We begin with a model from the series constructed by Dejonghe (1987). These are parametrized by a single dimensionless parameter, $q$, in the range $(-\infty < q \le 2)$. Increasing positive values correspond to models with increasing radial anisotropy, while increasingly negative values correspond to those with increasingly tangential anisotropy. The Einstein–Plummer model, with all stars on circular orbits, corresponds to the limit $q \to -\infty$.

(ii) The velocities are adjusted so that a fraction $\alpha$ of the particles with negative azimuthal velocity $v_\phi$ (in spherical coordinates $r, \theta, \phi$) are changed from negative to positive in a certain region of phase space, as described below. Thus, $\alpha = 0$ corresponds to Dejonghe's model, and when $\alpha = 1$ all particles have positive azimuthal velocity. We refer to this change in the velocities as the action of 'Lynden-Bell's demon' (Lynden-Bell 1962), and we use the terms 'prograde' and 'retrograde' to mean that $v_\phi > 0$ or $<0$, respectively. (The demon can be applied in reverse, as in (b) below.)

(a) **Basic models:** In the simplest series of models, the demon is allowed to work only on a certain fraction of the mass ($0 \le \gamma_1 \le 1$) when the particles are ordered by energy. Thus, $\gamma_1 = 0.2$ means that the Lynden-Bell demon works on the 20 per cent most bound particles; $\gamma_1 = 0$ gives a Plummer–Dejonghe model, while when $\gamma_1 = 1$ and $\alpha = 1$, all particles are given prograde rotation. These models are intended to give insight into the range of radii that are active in the instabilities observed in them.

(b) **High-shear models:** In the next series of models, the above procedure is applied to all particles (or rather to a fraction $\alpha$ of them) in the most bound fraction $\gamma_2$, but in the rest, i.e. the least bound fraction $1 - \gamma_2$ of them, the Lynden-Bell demon acts in reverse, reversing $v_\phi$ for all prograde stars (or a fraction $\alpha$ of them). These models are intended to maximize the shear in the rotation profile (Section 2.3), as this may be a factor in the existence of instabilities.

(c) **Low-inclination models:** Third, we consider a series of models in which the demon works only on stars of low inclination $i < \gamma_3\pi/2$, where $i$ is defined in the sense of the inclination of the orbital plane to the equatorial plane $z = 0$, so that $0 \le i \le \pi/2$. A motivation here is to investigate to what extent an instability can be thought of as being driven by motions in a thick conical disc.

(d) **High-inclination models:** Lastly, we have a set of models in which the Lynden-Bell demon works only on stars of high inclination $i > \gamma_4\pi/2$. Such models may give additional indications of the geometrical and kinematic characteristics of orbits that contribute strongly to instability.

These are the only initial models considered in this paper, but many variants are possible. For example, in the high-shear models, the axis of rotation of the least bound stars may be chosen differently from that of the most bound stars. Many such models can be explored by setting up initial conditions with the software PlummerPlus.<fnref n="1" />

<footnote n="1">¹https://github.com/pgbreen/PlummerPlus</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Rotation curves of six ‘basic’ models with $q = -6$, $\alpha = 1$ and values of $\gamma_1 = 0(0.2)1$. Here, $\gamma_1$ is the fraction of stars, ordered by energy, below which the Lynden-Bell demon is applied. The ordinate is the mean azimuthal velocity in a cylindrical shell of radius $R$ (the abscissa). Each of the 200 shells contains 500 particles. The six plots correspond to application of the Lynden-Bell demon (with various values of $\gamma_1$) to the same realization of a Plummer–Dejonghe model; therefore the fluctuations in successive models are correlated.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Rotation curves of six ‘high-shear’ models with $q = -6$, $\alpha = 1$ and values of $\gamma_2 = 0(0.2)1$. As Fig. 1, except that here $\gamma_2$ is the fraction of stars, ordered by energy, below which the Lynden-Bell demon is applied in the sense of prograde rotation, while above that energy it is applied in such a way as to give retrograde rotation.</figcaption></figure>

## 2.3 Characterization of rotation

There are several ways in which the rotation of these models can be measured. The first is a rotation profile, which can be defined in different ways, e.g. in spherical or cylindrical geometry, of which we concentrate on the latter. It is easily calculated from any model by first calculating the average rotation speed $\bar{v}_{\phi}(R) = \langle v_{\phi} \rangle$ as the average azimuthal velocity in a cylindrical shell of cylindrical radius $R$, and a thickness chosen as a compromise between sufficient spatial resolution and sufficient statistical reliability. From an $N$-body simulation, the corresponding quantity is easy to compute on a grid of $R$. Several examples are given below (Section 2.4).

For a global measure of rotation, we compute the rotational kinetic energy, but for this purpose define the average rotation speed on a circular ring of given $R$ and $z$ and, again, suitably chosen bin size. Then the rotational kinetic energy is defined to be

$$T = \frac{1}{2} \int \rho(R, z) \bar{v}_{\phi}^2(R, z) 2\pi R dR dz, \quad (1)$$

the integral being evaluated in practice by summing over the chosen grid. For our purpose, a fixed grid in $R$ with 100 bins, and a grid in $z$ adapted to contain approximately equal numbers of particles, was adopted. Of particular importance in this field, because of the pioneering work of Ostriker & Peebles (1973), is the dimensionless parameter $T/W$, where $W$ is the binding energy of the system.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Rotation curves of six ‘low-inclination’ models with $q = -6$, $\alpha = 1$ and values of $\gamma_3 = 0(0.2)1$. As described in the text, the Lynden-Bell demon is applied only to low-inclination orbits, with $i < \gamma_3\pi/2$. The ordinate is the mean azimuthal velocity in a cylindrical shell of radius $R$ (the abscissa). Each of the 200 shells contains 500 particles.</figcaption></figure>

## 2.4 Rotation profiles

Fig. 1 illustrates various rotation curves for the set of ‘basic’ models (Section 2.2). In this figure, the uppermost (blue) curve shows the rotation curve of a model where the Lynden-Bell demon is applied to all stars, whereas the lowest shows the rotation curve of the original Plummer model (anisotropy in all the models being modified by setting $q = -6$), and this is non-zero only because the profile has been determined from an $N$-body model. The green curve, corresponding to $\gamma_1 = 0.2$, shows the result when the demon is applied only to the

20 per cent most bound particles. Other choices of $\alpha$ simply alter the ordinate in proportion to $\alpha$.

Fig. 2 gives similar plots for the ‘high-shear’ (counter-rotating) models, in which $\gamma_2$ is the fraction of stars, ordered by energy, below which 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 differences between the maximum and minimum (largest retrograde) velocities are about twice as large.

Next, we can consider the effect of the inclination of the rotating orbits, by applying the Lynden-Bell demon only to orbits of low inclination; specifically, to stars such that $i < \gamma_3\pi/2$. The resulting rotation curves are illustrated in Fig. 3. At small cylindrical radius, this procedure greatly reduces the amount of rotation (unless $\gamma_3$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Rotation curves of six ‘high-inclination’ models with $q = -6$, $\alpha = 1$ and values of $\gamma_4 = 0(0.2)1$. The Lynden-Bell demon is applied only to high-inclination orbits, with $i > \gamma_4\pi/2$. The ordinate is the mean azimuthal velocity in a cylindrical shell of radius $R$ (the abscissa). Each shell contains 500 particles.</figcaption></figure>

is very large), since most stars at small cylindrical radius must lie on high-inclination orbits. The effect at large radius is much smaller, unless the critical value $\gamma_3$ is very small. In a complementary fashion, we can also choose to apply the demon only to high-inclination orbits, such that $i > \gamma_4\pi/2$, which leads to the rotation curves in Fig. 4.

## 2.5 Analysis of results

Our method of searching for a bar, in a spherical model which is rotating about the z-axis, is described in Rozier et al. (2019), along with studies of its accuracy.<fnref n="2" /> Briefly, a bar mode in the x, y plane is measured by computing the following integral, which uses cylindrical polar coordinates:

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

where $m = 2$ for a bar, and the integral is replaced by a sum of delta functions, one for each particle. Then the growth rate is computed by fitting an exponential<fnref n="3" /> to $|C_m(t)|$ over a suitable time interval. Actually, this is a relatively insensitive measure for non-rotating systems ($\alpha = 0$ in this work), as the bar (if there is one) need not form in the x, y plane. Therefore, non-rotating systems ($\alpha = 0$) can be measured by computing the corresponding integrals involving spherical harmonics of the order of $l = 2$, but this is not what is done in this study.

# 3 RESULTS

## 3.1 Basic models (prograde rotation)

We begin (Section 2.2) with the models in which prograde rotation is added at energies less than (more bound than) a certain fraction $\gamma_1$ of mass. Let us consider first the fiducial case $q = -6$. Fig. 5 shows the growth rate $\eta$ as a function of the Ostriker-Peebles parameter $T/W$, where $T$ is the rotational (mean flow) part of the kinetic energy.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Growth rate versus Ostriker-Peebles parameter for six sets of models with $q = -6$, $0 \le \alpha \le 1$ and various $\gamma_1$. Note that two realizations have been computed for the case $\gamma_1 = 1$, which gives an impression of the accuracy with which the growth rate has been estimated.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Growth rate for the maximally rotating model at each $\gamma_1$, for several values of $q$. We have avoided the isotropic case $q = 0$, as it is known (Rozier et al. 2019) that any instabilities here could belong to a different family from those at smaller (more negative) values of $q$. The result of two independent simulations is given for $\gamma_1 = 1$.</figcaption></figure>

Note that results for $\gamma_1 = 0$ are omitted, as these models are stable (in the sense that $\eta$ is not significantly non-zero), according to Rozier et al. (2019). Results for $\gamma_1 = 1$ are duplicated, i.e. two independent realizations are shown for each value of $\alpha$, to give an indication of the numerical uncertainty in the measured growth rate. Each sequence of symbols corresponds to a given value of $\gamma_1$, and one moves to the right within a given sequence as $\alpha$ increases. This increase in $\eta$ is satisfactorily fitted by a cubic power law in $T/W$, or equivalently

$$\eta(q, \gamma_1, \alpha) = \eta_{\max}(q, \gamma_1)\alpha^6, \quad (3)$$

a remarkably steep dependence.

Because of the scaling with $T/W$, we shall now only consider the function $\eta_{\max}(q, \gamma_1)$, results for which are shown in Fig. 6. Results of two models are given at $\gamma_1 = 1$, again to give an impression of the reproducibility of the results for independent initial conditions. Even the large fluctuations (such as the maximum at $\gamma_1 = 0.8$ in the lowest

<footnote n="2">²See also the caption to Fig. 5.</footnote>

<footnote n="3">³A logistic fit is used in Appendix A.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Growth rate for the maximally rotating model at each $\gamma_1$, for several values of $q$, plotted against the value of the rotation speed at radius $R = 0.5$.</figcaption></figure>

curve) may be explicable in terms of the choice of initial conditions, but trends in the data cannot be so discounted. The qualitative point that one may take from this plot is the fact that the growth rate does not keep increasing as the fraction $\gamma_1$ increases. Rather, the growth rate saturates when roughly the inner half of the system (ordered by energy) has been set to rotate in one sense.<fnref n="4" /> Turning back to Fig. 1, in particular the rotation curves for $\gamma_1 = 0.6, 0.8$, and 1, the large differences in these curves, which are particularly evident outside $r \sim 0.5$, make little difference to the growth rate. By contrast, the rotation curves for $\gamma_1 = 0.2$ and 0.4 change dramatically inside $r \sim 0.5$, as do the growth rates (Fig. 6). This region ($r \lesssim 0.5$), then, appears to be the region within which one must seek the mechanism for the instability, and any quantitative, phenomenological model for it.

To make this discussion more quantitative, in Fig. 7 we show how the maximum growth rate (as plotted in Fig. 6) varies with the value $\bar{v}_{\phi}(0.5)$. This plot shows that this value of $\bar{v}_{\phi}$ is a roughly linear predictor of the growth rate, though for very negative values of $q$ it does not predict well the saturation of the growth rate at the largest values of $\bar{v}_{\phi}(0.5)$. 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 $\alpha = 1$, i.e. the Lynden-Bell demon acts at its maximum level. But now the entire data set (all values of $q \le 0$, $\alpha$ and $\gamma_1$; see Fig. 8) can be grasped in outline by means of the foregoing analysis. The axes correspond to the Ostriker-Peebles stability parameter and the azimuthal streaming speed at the virial radius, and the size of the symbol encodes the growth rate of the bar. The region of the diagram containing the fastest growth rate is particularly interesting: although the classical Ostriker-Peebles criterion (for flattened systems) $T/W \gtrsim 0.14$ corresponds rather well with the boundary of this region when $\bar{v}_{\phi}(1)$ is large, we now find strongly unstable models when $T/W$ is much smaller, provided that $\bar{v}_{\phi}(1)$ is small. These are models in which $\gamma_1$ is sufficiently small that the azimuthal speed has dropped almost to zero by the virial radius (cf. Fig. 1). In fact the evident, roughly horizontal groupings in this

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Scatterplot of all results on models with prograde rotation, in the plane of $T/W$, $\bar{v}_{\phi}(1)$, i.e. the Ostriker-Peebles parameter and the azimuthal speed at the virial radius. The growth rate is coded in the size of the plotted points; the smallest symbols indicate growth rates less than 0.05.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Growth rate for the maximally rotating (and counter-rotating) model at each $\gamma_2$, for several values of $q$. For the reason stated in the text, the results for $\gamma_2 = 1$ have been duplicated for $\gamma_2 = 0$ (though new realizations would have led to slightly different results).</figcaption></figure>

diagram correspond to the two or three highest values of $\gamma_1$ included in our numerical survey.

## 3.2 High-shear models (mixed pro- and retrograde rotation)

For this sequence of models, we assume that the scaling with $T/W$ is similar to that in the purely prograde case, and proceed directly to the results for $\alpha = 1$. The growth rate is shown in Fig. 9. Qualitatively the results behave quite similarly to those for the purely prograde models, and even the maximum growth rates (which occur for models with smaller, i.e. more negative, values of $q$, and then at intermediate values of $\gamma_2 \simeq 0.5$) are not as dissimilar as the more extreme gradient in the velocity curve (compare Figs 1 and 2) might lead one to expect. Nevertheless one notices even more strongly that the maximum growth rate occurs at intermediate values of $\gamma_2$ (cf. footnote 4); furthermore the growth rate there is even larger than in Fig. 6. These remarks suggest that the maximum growth rate depends

<footnote n="4">⁴Actually there is some evidence from this figure that the saturation level is a little lower than the peak level, as one notices particularly for the most negative values of $q$.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Scatterplot of all results on models with pro- and retrograde rotation, in the plane of $T/W$, $|\bar{v}_{\phi}(1)|$, i.e. the Ostriker-Peebles parameter and the absolute value of azimuthal speed at the virial radius. The growth rate is coded in the size of the plotted points; the smallest symbols indicate growth rates less than 0.05.</figcaption></figure>

on the rotation of the outer (least bound) parts of the system: as $\gamma_1$ increases beyond 0.5, these parts start to rotate, and the growth rate is depressed slightly, while as $\gamma_2$ increases beyond 0.5 (Fig. 9) the outer regions switch from retro- to prograde rotation, and the drop in growth rate is more pronounced and starts from a higher maximum.

In Fig. 9, there are significant changes in the rate of growth at the left, around $\gamma_2 = 0.2$, where the values of $\eta_{\text{max}}$ have a very pronounced minimum, increasing again as $\gamma_2$ decreases further. In fact, for $\gamma_2 = 0$ the growth rates must be the same (within the error of the estimates) as for $\gamma_2 = 1$, as these two cases have the same rotation curve except for the sense of rotation.

If we contrast the great difference in $\eta_{\text{max}}$ between $\gamma_2 = 0$ and $\gamma_2 = 0.2$ with the small change in the rotation velocity at $R = 1$, it is clear that $\bar{v}_{\phi}(1)$ (coupled with $T/W$) cannot be as simple a predictor of the growth rate as in the case of the prograde models (Fig. 8). In Fig. 10, we show all results for the counter-rotating models, except that now the ordinate is the absolute value of $\bar{v}_{\phi}(1)$. The result is qualitatively similar to that for the prograde models, but more noticeable here is the fact that there are some very rapidly rotating models (high $T/W$ and $|\bar{v}_{\phi}(1)|$) with relatively low growth rates. We now investigate the rotation curves of such models to try to identify distinctive features of their rotation curves.

In fact, the models with low growth rate, high $T/W$ and high $|\bar{v}_{\phi}(1)|$ tend to have small (very negative) $q$ and $\gamma_2 \sim 0.2$. Such models are also very distinctive in Fig. 9, as discussed above. What distinguishes the rotation curve for these models, as can be seen in Fig. 2, is that the interior radii of a $\gamma_2 = 0.2$ model are hardly rotating, at least by comparison with radii above about 0.5. We can conjecture that the growth rate is small when the rotation speed $\bar{v}_{\phi}(0.2)$ is small. Therefore, to separate such models better than in Fig. 10, we simply multiply the abscissa by $|\bar{v}_{\phi}(0.2)|$. The result is shown in Fig. 11. There is now a much better separation of slowly and rapidly growing modes, and this choice of parameters has comparable success in the basic models of Section 3.1. These results confirm our finding (in that subsection) that the kinematics of the region inside $R \simeq 0.5$ is of particular importance for the presence of instability.

While this discussion has been focused on the growth rate, it is also interesting to measure the pattern speed in relation to the coexistence

<figure><img src="image_14.png" /><figcaption>**Figure 11.** As Fig. 10, except that the abscissa is multiplied by the mean azimuthal speed at radius $R = 0.2$, to separate models that are rotating relatively slowly at small radii. The growth rate is coded in the size of the plotted points; the smallest symbols indicate growth rates less than 0.05.</figcaption></figure>

of pro- and counter-rotation at different radii. The model with $\gamma_2 = 0$ is completely counter-rotating (Fig. 2), and for $q = -16$ (which we focus on in this paragraph) the pattern speed is also negative: $\Omega_b \simeq -0.7$.<fnref n="5" /> The pattern speed reduces in magnitude to $-0.5$ at $\gamma_2 = 0.1$, but is about +0.7 (i.e. prograde) throughout the range $0.3 \le \gamma_2 \le 1$. At $\gamma_2 = 0.2$, the instability is complicated, and no pattern speed could be measured confidently: up until $t \simeq 40$ a weak bar grows with a fast pattern speed around $\Omega_b = -1.7$, but this gives way by $t \simeq 60$ to a bar with faster growth and $\Omega_b \simeq -0.4$.

### 3.3 Models with prograde rotation at low or high inclination

More briefly, we consider the maximal growth rate (i.e. for $\alpha = 1$) for models in which rotation is given only to orbits of low inclination, i.e. those with rotation curves exemplified by Fig. 3, or only to orbits of high inclination (Fig. 4). The results are shown in Fig. 12, with the low-inclination series on the right. As $\gamma_3$ decreases from 1 to about 0.8, the growth rate drops by about a factor of 2, showing that stars with inclination above about 70° are very important in the instability. (This does not mean that stars of low inclination are irrelevant, but they are not sufficient of themselves for rapidly growing instability.) A complementary conclusion is drawn from the high-inclination models, whose growth rate drops by a comparable factor if $\gamma_4$ increases from 0 to about 0.3, i.e. if one excludes orbits with inclination less than about 30°. Thus, the presence of both high- and low-inclination stars is necessary for rapid growth, but neither group by itself is sufficient.

### 3.4 Spatial-temporal structure of the instability

#### 3.4.1 Spatial structure

In this subsection, we focus on the basic case $q = -16$, $\alpha = \gamma_1 = 1$, where the effects of rotation are strong, and apparently representative of all the very tangentially anisotropic models we have studied,

<footnote n="5">⁵Measurements reported in this paragraph were taken in a time interval around $7 < t < 21$, from the time of the clear emergence of the bar up to a time which is about half of the time when it reaches its maximum amplitude. The values are thought to be accurate to better than 10 per cent.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 12.** As Fig. 9, except that the Lynden-Bell demon is applied at only low or high inclinations, as determined by the parameters $\gamma_3$ and $\gamma_4$, respectively.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Projected density, in the plane $z=0$ of the bar in an $N$-body model with $N = 16384$ at $t = 20$HU. The result has been obtained by subtracting the projected density at $t=0$, and therefore shows only the perturbation. The data were gridded with a pixel size $0.2^2$HU, and contours were smoothed by interpolating to a grid 9 times finer. The scale gives the number per pixel.</figcaption></figure>

including the Einstein sphere ($q = -\infty$). Our aim is to present briefly some information on the spatial extent and shape of the bar. For other values of $q$, representations of the density are shown in Rozier et al. (2019, figs 5 and 6 for $q = -6$, and fig. 10 for other $q < 0$).

In $N$-body simulations with $N = 16384$, which were conducted using `NBODY6`, the bar is detectable (by Fourier analysis) only inside the cylinder $R \lesssim 1.5$, $|z| \lesssim 0.9$. Fig. 13 displays the projected density in the plane $z=0$ at time $t = 20$ (in this model the bar reaches its maximum amplitude at about $t = 40$). This can be compared with the result in Rozier et al. (2019, the uppermost inset in their fig. 10), obtained by the matrix method (Section 4.2.1). The maximum projected radius in Fig. 13 (at the corners) is nearly the same as the maximum cylindrical radius quoted above, and faint extensions can be seen at these radii. Also of interest here are the parts at right angles to the bar with negative projected density. In the $N$-body model at time $t = 20$, the core radius has decreased from an initial value of about 0.32HU to a value of about 0.26HU, presumably because of

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Growth rate and pattern speed of direct $N$-body models with $N = 32768$, $q = -16$, $\gamma_1 = 1$ for various values of $\alpha$. For $\alpha = 0, 0.2$ results are obtained by combining 10 models each, for $\alpha = 0.4$ there were eight, and there was one model each for the three largest values of $\alpha$. For $\alpha = 0$, the confidence interval on the growth rate extends to negative values, and only the upper limit is shown. Also shown, to illustrate the theoretical discussion of Section 5.1, is the maximum of $\Omega_2 - \Omega_1/2$ for nearly circular orbits, which is related to the inner Lindblad resonance (ILR).</figcaption></figure>

early core collapse. Possibly for this reason, the projected density at the centre of the bar is above zero.

### 3.4.2 Growth rate and pattern speed

The results that follow are based on a dedicated set of direct $N$-body models with $N = 32768$ particles, with anisotropy parameter $q = -16$. They are basic models of the type discussed in Section 3.1 with $\gamma_1 = 1$, i.e. the Lynden-Bell demon is applied everywhere. Our aim is to discuss the dependence of the time scales of the bar instability on the rotation parameter $\alpha$, with a view to informing the discussion of Section 5.1.2. Detailed discussion of the methods of analysing the runs are given in Appendix A. Results are given in Fig. 14.

While details of the analysis and interpretation of this plot are given in later sections and the appendix, explanation is required here for the fact that any pattern speed at all has been measured for the non-rotating case $\alpha = 0$. We emphasize that no growth could be measured with confidence; all that is shown in Fig. 14 is an upper limit. But the analysis of the $m = 2$ Fourier coefficient $C_2$ (see equation 2), using the autocorrelation of $\arg(C_2)$, clearly demonstrated the presence of oscillation with the stated pattern speed. The figure supports the view that the bar instability in models with $\alpha > 0$ connects continuously with the dynamical behaviour of the non-rotating model, but does not imply that the non-rotating model is unstable.

## 4 RESONANCES

### 4.1 Resonances in tangentially biased spherical systems

The role of resonances in the occurrence of bar instabilities in discs is sufficiently well known to be textbook material (e.g. Binney & Tremaine 2008). Less familiar, however, is their role in spherical systems, despite the substantial literature referred to in Sections 1 and 2.1. In this introductory subsection, we focus on (very) tangentially anisotropic systems, where the radial motions are small, and we can
<table><caption>**Table 1.** Resonances in the leading-order perturbations of epicyclic motions in a rotating spherical model. Names that are not in common use are given in slanted type. Here, the orbital inclination *i* is used in its conventional sense, with retrograde orbits having $\pi/2 < i \le \pi$.</caption><thead><tr><th>Ω<sub>b</sub> =</th><th>i-dependence</th><th>Name</th><th>(n<sub>1</sub>, n<sub>2</sub>)</th></tr></thead><tbody><tr><td>0</td><td>sin<sup>2</sup>i</td><td>tumbling</td><td>(0,0)</td></tr><tr><td>Ω</td><td>(1 + cos i)<sup>2</sup></td><td>corotation</td><td>(0,2)</td></tr><tr><td>Ω − κ/2</td><td>(1 + cos i)<sup>2</sup></td><td>ILR</td><td>(−1, 2)</td></tr><tr><td>Ω + κ/2</td><td>(1 + cos i)<sup>2</sup></td><td>OLR</td><td>(1,2)</td></tr><tr><td>±κ/2</td><td>sin<sup>2</sup>i</td><td>epicyclic</td><td>(±1, 0)</td></tr><tr><td>−Ω</td><td>(1 − cos i)<sup>2</sup></td><td>reverse corotation</td><td>(0, −2)</td></tr><tr><td>−Ω + κ/2</td><td>(1 − cos i)<sup>2</sup></td><td>reverse ILR</td><td>(1, −2)</td></tr><tr><td>−Ω − κ/2</td><td>(1 − cos i)<sup>2</sup></td><td>reverse OLR</td><td>(−1, −2)</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Resonances related to circular motions in a Plummer potential perturbed by a bar. Also shown is the approximate pattern speed of a bar in a rotating model with $q = -16$, $\alpha = 1$ and $0.4$.</figcaption></figure>

employ the epicyclic approximation. Then the main resonances can be identified by a quite elementary calculation (Appendix B); it leads to the results given in Table 1, which we now explain.

We suppose there is a bar perturbation, whose pattern rotates with angular speed $\Omega_b$. Then the familiar corotation resonance (line 2) is given by $\Omega_b = \Omega$, where $\Omega$ is the circular angular velocity. Similarly, the Inner and Outer Lindblad Resonances (ILR, line 3; OLR, line 4) are defined by $\Omega_b = \Omega \mp \kappa/2$, where $\kappa$ is the epicyclic frequency. What is more unfamiliar is their inclination-dependence; and what we have termed *reverse* resonances would be the normal corotation and Lindblad resonances for a bar with retrograde pattern speed (as one can see from the *i*-dependence). Finally, there are what we have termed *epicyclic* and *tumbling* resonances, which are strongest at high inclinations (near $\pi/2$).

In a more general setting (Section 4.2), these resonance relations are written as $2\Omega_b = n_1\Omega_1 + n_2\Omega_2$, where $n_1, n_2$ are integers, $\Omega_1$ is the frequency of radial motions, and $\Omega_2$ is the mean longitudinal frequency, which generalize $\kappa$ and $\Omega$, respectively. These number pairs are included in the last column of the table.

The simplest situation in which to visualize the possible location of these resonances is the limiting case $q \to -\infty$ (the 'Einstein sphere'), when the frequencies $\Omega_1, \Omega_2$ become $\kappa, \Omega$. Combinations of these frequencies are plotted in Fig. 15. Also shown is the observed pattern speed (Section 3.4.2) found for the model with $q = -16$ in

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Distribution of frequencies $\Omega_1, \Omega_2$ at time 15 in an N-body model with $N = 32768$, $q = -16$, $\alpha = 1$, $\gamma_1 = 1$, along with the loci of five resonances. At this time the model exhibits a growing bar (Section 4.3), and the pattern speed was determined in its early growth phase. OLR/ILR: Outer/Inner Lindblad Resonance, CR: Corotation Resonance, ER: epicyclic resonance (Table 1), rILR: reverse ILR.</figcaption></figure>

the cases $\alpha = 0.4, 1$. Here, the anisotropy is sufficiently tangentially biased that the Einstein sphere is approximately representative.

In a later section (Section 4.3), we shall present some evidence of the influence of resonances on the behaviour of particles in N-body models. Here, for one model, we simply present in Fig. 16 the distribution of particles in the space of the frequencies $\Omega_1, \Omega_2$, where their proximity to the various resonances can be appreciated visually. For the initial conditions of a model with $N = 32768$ particles, $q = -16$ and $\alpha = 1$ ($\gamma_1 = 1$ in Section 2.2), the values of $\Omega_{1,2}$ were calculated by integrating the equations of motion in the Plummer potential.

The distribution is narrow because the anisotropy is quite extreme ($q = -16$) (in the limit $q \to -\infty$ the distribution would have been 1D). The ILR and its retrograde counterpart appear to involve no particles, but we shall see (Section 4.2.2) that the former still appears to play a central role in the dynamics, and present below a possible reason for this. There are, however, particles lying close to the other three resonances, and empirical evidence for their influence is studied in Section 4.3. We can see immediately, however, that these particles are confined to the low-frequency part of the distribution, and this means that they affect stars at relatively large radius. With the exception of the 'tumbling' resonance, the remaining resonances in Table 1 appear to play no role.

In this section, we identified resonances by examining the Fourier decomposition of the perturbing potential. In first-order orbital perturbation theory, a given Fourier term gives rise to perturbations containing a resonant denominator $\omega - n_1\Omega_1 - n_2\Omega_2$, where $\omega$ is the frequency of the explicit time-dependence of the perturbing potential. If resonance is strictly defined as the vanishing of this denominator (Jalali & Hunter 2005, p.810) then resonances do not occur if the perturbation grows, but 'only near-resonances at which the denominators [...] are small' (if the growth-rate is small). But these authors also (their p.819) point to a situation, closely analogous to ours, in which the Fourier component corresponding to the ILR is important even though no particle is in strict or near-resonance. It is a situation in which 'the denominator [...] is never large for any orbit', an apt description for the fact that the line labelled 'ILR' in Fig. 16 is approximately parallel to the distribution of points. Such a Fourier
term, it seems, may have a comparable importance to a resonant term in which a small fraction of particles are in 'near-resonance' but for most particles the denominator is large. Though it can be argued that there is no resonance at all if the denominator is not small compared with typical orbital frequencies, we prefer to regard the ILR as a weak resonance affecting most, if not all particles. In general, in what follows we take a liberal view of what is meant by 'resonance'.

## 4.2 Analysing the role of resonances with the matrix method

Not all theoretical discussions of bar instabilities in discs seem readily adaptable to spherical systems, and here we follow specifically the arguments proposed by Palmer et al. (1989) and Allen et al. (1992). These two papers were concerned largely with two instabilities. In the earlier paper, they developed a theory of what they termed the circular orbit instability, generated by the ILR, and in the later paper they turned to a so-called tumbling instability, generated by what we have (for that reason) termed the tumbling resonance (Table 1). Both papers proceed from the study of the linearized, orbit-averaged collisionless Boltzmann-Poisson system. More precisely, they demonstrate the existence of neutral (or unstable) solutions to the resulting eigensystem by only keeping terms that correspond to a particular resonance (or a particular set of resonances) between the instability and the orbital frequencies in the system. We propose to use a similar method, and try to test whether these mechanisms can be responsible for instabilities in the systems studied in this paper.

### 4.2.1 Method

In order to investigate in detail the role of resonances in the development of instabilities in tangentially anisotropic, rotating systems, we extended the version of the matrix method formulated in Rozier et al. (2019) to the case of rotating spheres. The matrix method identifies unstable modes as eigenvectors (with an eigenvalue equal to 1) of the response matrix

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

where the integral runs over the action space ($J_r, L, L_z$) consisting of the radial action, the total angular momentum, and its z-component, respectively; $F$ is the phase-space distribution function of the system; $\Omega = (\Omega_r, \Omega_\phi, 0)$ are the frequencies of a given orbit; $\psi_n^{(p)}$ is a Fourier-transformed element of a bi-orthogonal potential-density basis; $n$ is a three-vector of integers; and $\omega = \omega_0 + i\eta$ is the (complex) frequency of the instability (oscillation frequency and growth rate). Examining the integrand in equation (4), we can see that the denominator is resonant, in the sense that it can take very small values in some parts of action space. Those orbits where $\omega - n \cdot \mathbf{\Omega}(\mathbf{J})$ is close to 0 are expected to be decisive in the development of the instability, since they contribute the most in the response matrix. We then call the 3D vector of integers a resonance vector: at a given $n$, the corresponding integral will mainly involve the orbits that are at the $n^{th}$ resonance with the instability.

As a consequence, each integral term in the sum of equation (4) can be regarded as the contribution to the response matrix from the orbits that are at the $n^{th}$ resonance with the instability. So, in order to investigate the role of each resonance in the development of instabilities, one can separately study the influence of each corresponding term in the matrix. We decided to follow this principle in order to study the interplay between resonances and instabilities in tangentially biased, rotating systems. More specifically, instead of

just comparing the amplitude of the resonance terms in equation (4), we opted for the full mode search developed in Rozier et al. (2019), applied to response matrices in which only a subset of the resonance terms were kept: if $\mathcal{I}$ is a sub-ensemble of $\mathbb{Z}^3$, we searched for modes assuming the response matrix is restricted to

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

Such an approach has already been used in a study of disc instabilities by Polyachenko (2005).

An important constraint on the ensemble of resonances to consider comes from the fact that $\hat{M}_{pq}(\omega) \propto \delta_{mp}^{n_3} \delta_{mq}^{n_3}$. Since we are focusing on 2-fold symmetric instabilities,<fnref n="6" /> we only compute matrix terms with $m^p = m^q = 2$, hence the sum in equation (5) is restricted to $n_3 = 2$. Consequently, in the following paragraphs, we will consider $\mathcal{I}$ as a sub-ensemble of $\mathbb{Z}^2$, and define $\tilde{n} = (n_1, n_2)$.

Using this notation, we can label some crucial resonances by their resonance vector $\tilde{n}: (-1,2)$ corresponds to the ILR, (0,2) is the corotation resonance, (1,2) is the OLR, and so on (Table 1). Each resonance should be understood through the relative value of the vector elements, so for each resonance vector $\tilde{n}$, its opposite $-\tilde{n}$ corresponds to the same kind of resonance.<fnref n="7" /> Hence, whenever a resonance vector will be included in the ensemble $\mathcal{I}$, it will be assumed that its opposite is also included. As identified in Rozier et al. (2019), there is little influence from $|n_1| > 2$ to the instabilities in tangentially anisotropic, rotating systems, so that we will restrict the study to the pairs $(n_1, n_2)$ with $|n_1| \le 2$. As a consequence, the case of the complete matrix has eight different resonance vectors: $\mathcal{I}_0 = \{(0, 0), (0, 2), (-1, 2), (1, 0), (1, 2), (-2, 2), (2, 0), (2, 2)\}$, the first five of which are listed in Table 1.

Finally, we draw attention to the fact that a major difference between rotating and non-rotating systems emerges in the case $\tilde{n} = (0, 0)$. The corresponding term in the response matrix (equation 4) is directly related to the $L_z$-gradient of the distribution function, and is independent of the orbital frequencies. We will refer to it as the tumbling term: while it does not correspond to any resonance involving orbital frequencies, it is shown in Allen et al. (1992) that this term is associated with the tumbling of orbital planes under the effect of a perturbation (see also Palmer 1994).

In summary, our aim is to investigate the role of resonances by examining the effect on the growth rate and pattern speed of removing some of the corresponding terms from the calculation. In order to strengthen the case for this method, we present in Appendix C the result of its application to the well-studied radial orbit instability regime found in Rozier et al. (2019). We now apply it to the tangentially biased rotating systems that are the focus of this paper, and concentrate on the models with $q = -16$ and $\alpha = 0.4, 1$.

### 4.2.2 Application to instability in tangentially biased systems

In order to identify which resonances are the most critical to the instabilities, we compared the results of the full matrix with results when a single resonance is removed. The importance of the removed resonance in the instability is then estimated through the change in

<footnote n="6">⁶Note that in principle, other instabilities could exist (i.e. with $m = 0$, giving axisymmetric oblate/prolate oscillations), however, they are beyond the scope of this paper.</footnote>

<footnote n="7">⁷Note that the role of opposite resonance vectors is not the same; they only correspond to the same kind of resonances. This is illustrated by the different inclination-dependence of resonances and their 'reverse' analogues in Table 1.</footnote>
<table><caption>**Table 2.** Values of the growth rate and pattern speed found through the restricted response matrix method applied to the $(q, \alpha) = (-16, 0.4)$ and $(q, \alpha) = (-16, 1)$ tangentially biased systems: comparison between the matrix of reference (complete) and the series of matrices obtained by removing a single resonance.</caption><thead><tr><th>$\mathcal{I}_0 \setminus \mathcal{I}$ (removed)</th><th>$\eta(0.4)$</th><th>$\Omega_b(0.4)$</th><th>$\eta(1)$</th><th>$\Omega_b(1)$</th></tr></thead><tbody><tr><td>$\emptyset$ (reference)</td><td>0.012</td><td>0.36</td><td>0.24</td><td>0.66</td></tr><tr><td>(-1, 2)</td><td>0.0002</td><td>0.13</td><td>0.0066</td><td>0.33</td></tr><tr><td>(1,2)</td><td>0.0082</td><td>0.33</td><td>0.088</td><td>0.60</td></tr><tr><td>(0,2)</td><td>0.0040</td><td>0.30</td><td>0.022</td><td>0.49</td></tr><tr><td>(0,0)</td><td>0.0073</td><td>0.33</td><td>0.044</td><td>0.40</td></tr><tr><td>(1,0)</td><td>0.0044</td><td>0.36</td><td>0.20</td><td>0.66</td></tr><tr><td>(-2, 2)</td><td>0.010</td><td>0.36</td><td>0.24</td><td>0.66</td></tr><tr><td>(2,2)</td><td>0.010</td><td>0.33</td><td>0.20</td><td>0.66</td></tr><tr><td>(2,0)</td><td>0.011</td><td>0.36</td><td>0.24</td><td>0.66</td></tr></tbody></table>

growth rate and oscillation frequency of the instabilities found with and without the resonant term.

Table 2 shows the values of the growth rates and pattern speeds that were obtained by the restricted matrix method in the two cases $\alpha = 0.4$ and $\alpha = 1$, alternatively removing each of the resonant terms that are present in the matrix of reference. Notice that the values of the growth rate and pattern speed of reference do not exactly correspond with those of Rozier et al. (2019): for the sake of numerical efficiency, we opted for fewer radial basis functions, so the mode reconstruction is not as accurate. However, the present reference case identifies the same instability, and we ensure the self consistency of the present method by performing matrix calculations with constant parameters, in particular the nature and number of radial basis functions. Similarly, the discrepancy in the measured growth rate of the $\alpha = 0.4$ cluster between simulations (see Fig. 14) and matrix method is probably due to the relatively low number of basis functions used.

Let us first focus on the $\alpha = 0.4$ case. Some resonances have little influence both on the growth rate and on the pattern speed of the mode: removing either of $(-2, 2)$, $(2,2)$ or $(2,0)$ does not impact the frequencies by more than 20 per cent.<fnref n="8" /> Some resonances have a more significant impact on the instability growth rate: removing a term among $\{(1, 2), (0, 2), (0, 0), (1, 0)\}$ depletes the growth rate by at least 30 per cent, and up to 70 per cent. The pattern speed, however, is barely impacted by the removal of any of these terms. It appears that the contribution to the instability from the resonances in these groups can be interpreted as follows: none of them has a fundamental role in the formation of the instability, however, they all contribute more or less to it by increasing its growth rate. This analysis does not apply to the ILR term $(-1, 2)$: when it is removed, the remaining instability has a growth rate about 50 times smaller than the reference one, while its pattern speed is depleted by about a factor 3. This indicates that the instability identified in the absence of the ILR term is of a different nature from the one identified in all the other cases.

The maximally rotating, $\alpha = 1$ cluster seems to present a more subtle behaviour. One can still identify a group of resonances that have little significance to the instability: $\{(1, 0), (-2, 2), (2, 2), (2, 0)\}$. The group of resonances that influence the growth rate but not the pattern speed is now reduced to one element: the removal of the OLR

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Evolution of the distribution of $\Omega_2$ during the early growth of a bar in an $N$-body model with $N = 32768$ particles, $q = -16$ and $\alpha = 1$. The black bar shows the variation of $\Omega_b$ (which decreases slightly) during the same interval of time. The green bars show the ranges of the epicyclic (right) and Outer Lindblad (left) resonances. For those resonances, the breadth of frequencies is mainly due to the width of the distribution of frequencies in the initial conditions (Fig. 16).</figcaption></figure>

term results in the depletion of the growth rate by a factor of 3 and of the pattern speed by about 10 per cent. A new category emerges: the tumbling and corotation terms are now associated with decrease of a factor, respectively, about 5 and 10 in the growth rate, and of, respectively, about 40 per cent and 25 per cent in the pattern speed. Finally, the ILR term is still the one bearing the most importance in the instability, as its removal leads to a depletion of a factor about 40 in the growth rate, and about 2 in the pattern speed.

The most important process that leads to instability seems, in both $\alpha = 0.4$ and $\alpha = 1$ cases, to be borne by the ILR term in the matrix. This is indicative of a similarity between the instabilities in both cases. However, the behaviour of the other terms changes when rotation increases: the importance of the $(-2, 2)$, $(2,2)$, $(2,0)$ and especially $(1,0)$ terms is lower at higher rotation, while corotation, the OLR and the tumbling terms gain importance with rotation. The tumbling term seems to gain a particular role in setting the high pattern speed of the instability at high $\alpha$.

Before we pass on from the matrix method, it is worth mentioning that we have focused here on the most rapidly growing mode. For a broad range of values of large $\alpha$ and negative $q$ (in models where the Lynden-Bell demon is applied everywhere), it coexists with a second mode that has smaller growth rate and pattern speed (by factors of the order of 5 and 2, respectively).

## 4.3 Resonances in an $N$-body model

In this section, we return to $N$-body modelling to examine evidence for resonant behaviour. In particular, we discuss direct $N$-body simulations of basic models with $\gamma_1 = 1$ (i.e. the Lynden-Bell demon is applied everywhere), especially the case of extreme tangential anisotropy ($q = -16$) and maximal rotation ($\alpha = 1$). This exhibits a rapidly growing bar that reaches maximum amplitude at about $t = 30$, but its growth is clearly non-linear by $t = 15$ (when a plot of the amplitude against time has an inflection from concave upwards to concave downwards). During this shorter interval the pattern speed $\Omega_b$ 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 azimuthal frequency $\Omega_2$ during the period of roughly linear behaviour. In most cases, it is difficult to distinguish features from

<footnote n="8">⁸It is not hard to show that these terms are smaller than the others by a factor at least of the order of the epicyclic amplitude, which is small in these tangentially anisotropic models.</footnote>
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Change in the distribution of $\Omega_2$ between $t=0$ and $t=15$, given separately for particles of high ($i > 60^{\circ}$; green) and low inclination (red); the thick black line gives the sum, i.e. all inclinations. The error bars on the latter give the average and $1\sigma$ standard error of the mean over 10 simulations. For the first two graphs the error bars would be larger by about 40 per cent. The horizontal bars at the foot have the same meaning as in Fig. 17, but depend on the bar pattern speed, which was determined for only one of the simulations. The binwidth is half that of Fig. 17.</figcaption></figure>

sampling fluctuations, and we shall improve the statistics below, but in view of the large numbers of particles per bin, it is hard to dispute that there is a significant increase around $\Omega_2 = 0.65$ and an equally significant decrease around 1.05 (the mid-point of the histogram bar is quoted). It is tempting to link the first of these with corotation. Again tentatively, the second could be associated with the 'epicyclic' resonance: in Fig. 16 the corresponding line crosses the upper margin of the distribution of points, where they are concentrated, at a frequency just below $\Omega_2 = 1$. Other features at higher frequency, if regarded as significant, correspond to none of the resonances listed in Table 1, but there is another at a lower frequency, about $\Omega_2 = 0.4$, which could very well correspond to the OLR, as can be seen from Fig. 16.

It is worth examining the distribution of those stars that contribute most to the evolution of $f(\Omega_2)$. Fig. 18 shows two plots corresponding, respectively, to low and high inclination, and also the plot for all inclinations (black). The red curve, for instance, shows the change from $t=0$ to $t=15$ in the number of particles in a bin (centred at the indicated values of $\Omega_2$), counting only those of inclination $i < 60^{\circ}$. Results from 10 independent simulations have been combined, and the bin size is half that of Fig. 17. The sum, plotted in black, adds considerable significant detail to what can be inferred from Fig. 17. In the features at $\Omega_2 \simeq 0.4$ and $0.65$, the change is predominantly due to particles of low inclination, but the opposite is true of the feature at $\Omega_2 \simeq 1.05$. Bearing in mind the inclination-dependence noted in Table 1, these results lend weight to the interpretations as behaviour at the outer Lindblad, corotation and epicyclic resonances, respectively.

Despite what has just been said about the epicyclic resonance, actually the resonant frequency is somewhat too small (see the position of the rightmost horizontal bar in the figures), and in this respect a better interpretation is offered by consideration of the ILR, the importance of which has already been established in Section 4.2. It is defined by $\Omega_b = \Omega_{ILR} \equiv \Omega_2 - \Omega_1/2$, and it is found that, in the limit $q \to -\infty$, the expression $\Omega_{ILR}$ reaches its maximum value at a radius where $\Omega_2 = 1.09$. Its importance is emphasized again in Section 5.1. A third possible interpretation is that the feature near $\Omega_2 = 1$ (from about 0.8 to 1.1) is composite.

As with the inclination, one may study the radial distribution in the bins which we have associated with resonance, but the radius is essentially given by the frequency $\Omega_2$, which is a decreasing function of radius. Indeed the bins at $\Omega_2 = 0.4, 0.65$, and 1.05 correspond to approximate radial ranges $2 < r < 2.8, 1.2 < r < 1.6$, and $0.6 < r < 0.8$, respectively.

What has been learned from studying the evolution of the distribution of $\Omega_2$ can also be attempted with the distribution of $\Omega_1$, related to the epicyclic resonance, and $\Omega_1 + 2\Omega_2$, involving the OLR. But again because of the narrowness of the distribution of the two basic frequencies, one sees evidence of rapid evolution at essentially the same corresponding frequencies, i.e. those which may be associated with all three resonances, plus some evidence for unattributable evolution at still larger frequencies.

# 5 DISCUSSION, SUMMARY, AND CONCLUSIONS

## 5.1 Destabilizing mechanisms in tangentially anisotropic, rotating spheres

### 5.1.1 A Theoretical framework for bar formation

This paper has presented a diversity of empirical information on bar formation in a tangentially anisotropic, rotating spherical system. To tie everything together, it is helpful to set our findings within a theoretical framework, and the one which we now present, though not original, has proved most fruitful in the exploratory aspects of our research.

Palmer et al. (1989) and Allen et al. (1992) exhibit two particular mechanisms that can give rise to linear instability in, respectively, tangentially anisotropic and rotating spherical stellar systems: the circular orbit instability and the tumbling instability. Their approach is via a careful analysis of the linearized orbit-averaged Boltzmann-Poisson system, in which the first step is to treat one resonance in isolation, and then all others are treated as perturbations. For both mechanisms, it is shown in the first step that the system would exhibit a particular set of neutral modes (i.e. modes with vanishing growth rate) in the absence of the perturbing resonances, but that such a mode may be destabilized by the perturbing resonant interactions. While the tumbling process is linked to the precession of orbital planes under the influence of a rotating bar, the circular orbit instability emerges from the influence of a bar on the apsidal precession of the orbits, entailing changes in $\Omega_{ILR}$ (defined in equation 6).

In the case of the circular orbit process, the characteristics of the instability are the following:

(i) The ILR is at the source of the neutral mode, hence it is the fundamental cause of the instability.

(ii) The mode rotates at a pattern speed larger than (but close to) the maximum value of the ILR frequency, defined by

$$ \Omega_{ILR} = \Omega - \frac{1}{2}\kappa. \qquad (6) $$

Fig. 15 shows this frequency in the present clusters as a function of radius.

(iii) Destabilization is ensured by other resonances, at the oscillation frequency of the mode.

A similar list of characteristics can be given for the tumbling instability:

(i) The tumbling process [corresponding to the (0,0) resonance vector] is at the source of the neutral mode.
(ii) The pattern speed is nearly proportional to the rotation parameter $\alpha$. This comes from the fact that the tumbling term in the eigen equation describing the neutral mode (see equation 8.31 in Palmer 1994) is proportional to $\alpha/\omega$.

(iii) Destabilization is ensured by other resonances, at the frequency of the mode.

Before passing on to a discussion of empirical evidence on the mechanisms of bar formation, we mention here one factor that supports the framework of the above theory: the presence of a second mode (see the last paragraph of Section 4.2.2). The main problem of the cited authors' approach, in which only one resonance is included in the first instance, exhibits not one but a sequence of neutral modes, and there is no reason to suppose that only one of them can be destabilized by perturbing resonances.

### 5.1.2 The role of tumbling and resonances: empirical evidence

The series of experiments using the matrix method that we reported in Section 4.2 helps in associating the instabilities we exhibited through N-body experiments (mainly in Section 3) with the processes of tumbling and circular orbit instability. We begin mainly with the matrix results, and then pass on to an interpretation of the N-body data.

The $(q, \alpha) = (-16, 0.4)$ case is particularly eloquent: it is obvious (from Table 2) that the instability mostly depends on the ILR; the pattern speed of the mode lies just above the maximum of $\Omega_{ILR}$ (see Fig. 15); and the role of other resonances seems to be to accumulate the growth rate. The $(q, \alpha) = (-16, 1)$ case seems to follow a similar trend, yet less tightly matching the characteristics of the circular orbit instability: the ILR term still bears most of the instability; the pattern speed is twice as high as the maximum of $\Omega_{ILR}$; and yet some other terms play an important role, in particular the one associated with the tumbling process. In contrast to the conclusions of Allen et al. (1992), it seems that the instabilities in the present tangentially anisotropic, rotating systems cannot be attributed to the tumbling instability alone. Perhaps the high values of the pattern speed and the growth rate may diminish the accuracy of their perturbative approach. The results of Table 2 rather point to a mixing of the circular orbit and the tumbling processes, the ILR still being the most important resonance in the creation of the instability.

It is at first sight paradoxical to assert that the ILR is central to the observed instabilities, when Fig. 16 shows that no particles in the models lie close to resonance. However, the restricted matrix method (Table 2) shows that, even if the pattern speed is far away from any value $\Omega_{ILR}$ can take in the system (more than twice its maximum value when $\alpha = 1$), the ILR is fundamental to the development of an instability in the system. Furthermore, there are indications in the N-body models that the ILR is playing a vital role: the frequency measurements presented in Figs 17 and 18 show strong activity in particles for which $\Omega_2 = 1.05$, which corresponds closely to the value of $\Omega_2$ at which $\Omega_{ILR}$ takes its maximum value ($\Omega_{ILR} = 0.29$). Finally on this point, we recall Fig. 14, which illustrates the 'floor' that the pattern speed reaches as the rotation parameter $\alpha$ decreases; this floor is near the maximum of $\Omega_{ILR}$, as expected in the scenario outlined in Section 5.1.1.

The foregoing remark, that a resonance may have a vital role even in the absence of resonant particles, echoes a point made near the end of Section 4.1, and should be borne in mind in assessment of the role of the tumbling instability. (This is formally the condition $\Omega_b = 0$, and is never satisfied in our models.) Consider in particular the region $\Omega_2 > 1.1$ in Figs 17 and 18, which corresponds to relatively low radii. It is populated with orbits that do not resonate with the pattern at any

of the usual low order resonances, but there is evidence in Fig. 18 that the variations in this inner region mostly affect high-inclination orbits. Therefore, we interpret the activity here as being due to the tumbling process, because of its inclination-dependence (Table 1).

The results of Sections 3.1 and 3.2 show that much importance attaches to these inner orbits: from the high-shear models, we know that flipping the orbits at $R < 0.4$ is enough to nearly suppress the instability, while building a larger and larger rotating core of radius $0.4 < R < 0.75$ increases the growth rate of the instability in the inner system up to its maximum value. By contrast, the results of these two sections also show that what happens at larger radii ($R > 0.5$, say) hardly matters: this region may rotate, counter-rotate, or not rotate at all. The importance of the most central orbits is also supported by the 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 $R = 0.45$ in the present units, the density in the bar has dropped to 10 per cent of its maximum value. All the evidence supports the view that the tumbling process, while not dominant in the generation of instability, is an important component of the mechanism that creates it.

It should be noted that the foregoing discussion refers to the cylindrical radius $R$. For high-inclination orbits, such as those most affected by the tumbling and epicyclic resonances, the spherical radius is significantly larger. Still, Fig. 18 shows that there is also some activity in the low-inclination orbits (in the inner region on which we have been focusing).

Low-inclination orbits have a particularly important role in the regions of the OLR and corotation (around $\Omega_2 = 0.4$ and $0.7$, respectively); these are orbits at relatively large radius. The matrix results however show that the role of these resonances is of less importance than those of the tumbling or the ILR. This is consistent with the stability results of Section 3.3: both high and low-inclination orbits matter, but the high inclination ones are more important than the low-inclination ones.

A final factor that hints at the significance of the tumbling resonance at high $\alpha$ is the theoretical result (mentioned in Section 5.1.1) that the pattern speed of the relevant neutral mode is proportional to $\alpha$. This suggests an explanation of the rising pattern speed for $\alpha \gtrsim 0.5$ in Fig. 14; this could not continue for small $\alpha$, as then the damping effect of the ILR would come into play (Palmer et al. 1989, p.1292). It is interesting to note also, from Table 2 in the case $\alpha = 1$, that removal of the (0,0) term reduces $\Omega_b$ to a value quite similar to the reference result for the case $\alpha = 0.4$ which, we have argued, is constrained by the maximum of $\Omega_{ILR}$.

To sum up, the present results point to the following scenario for the formation of instabilities in tangentially anisotropic, rotating systems: in the inner region of the cluster, the tumbling of orbital planes (especially for high-inclination orbits) and the drift in apsidal precession far from ILR cooperate to the formation of an unstable bar. The growth rate of this pattern is enhanced by resonance with outer orbits at the epicyclic resonance and (most noticeably) the corotation resonance and the OLR, in particular with the low-inclination orbits. Thus, there is a kind of cooperation between all the resonances, which diminishes high-inclination orbits in the inner regions, and enhances those of low inclination at larger radii. Both effects create the flattened structure of the bar.

## 5.2 Towards a general stability criterion

In their classic study, Ostriker & Peebles (1973) stated that discs become 'approximately stable' when $T/|W|$ falls to a value of $0.14 \pm 0.02$. The wording shows just how difficult it is to establish
stability through numerical methods, and the same problem arises
in our work. Furthermore, the emphasis in their work was on discs
(with or without haloes), whereas ours is mainly (but not exclusively)
on spherical models that rotate globally. Nevertheless, we take the
Ostriker & Peebles criterion as our starting point in this subsection.

5.2.1 General rotating models studied in this paper

Figs 8, 10, and 11 summarize the results of our main numerical
surveys, placing them in the space of growth rate, *T/W* and particular
circular velocities in the cluster. Even if there seems to be a trend
towards more rapid instability with increasing *T/W*, an important
fraction of our parameter space does not satisfy the Ostriker-Peebles
criterion. The present results simply add to the number of published
deviations from this law (see e.g. Zang & Hohl 1978; Aoki et al.
1979; Berman & Mark 1979; Efstathiou et al. 1982; Evans & Read
1998; Athanassoula 2008), and yet are among the few that exhibit
such deviations in purely spherical clusters.

As a consequence, it appears that the linear stability of spherical
systems cannot be uniquely determined by the ratio between global
kinetic energy in rotating motion and global gravitational potential
energy. Nor would the global angular momentum be more useful: our
simulations of clusters with counter-rotating components (high-shear
models, Section 3.2) imply the presence of fast growing instabilities
in clusters with arbitrarily low global angular momentum (when the
angular momentum of the inner prograde and outer retrograde parts
of the cluster compensate).

An important reservation should be made at this point, con-
cerning the presence of instabilities with low growth rates. On
one hand, in terms of measurements by analytical means (i.e.
the response matrix in the present case), the behaviour of the
analytical indicators at low growth rate gets very intricate (see
e.g. Pichon & Cannon 1997; Merritt 1999; Rozier et al. 2019),
which makes it complicated to measure growth rates accurately.
On the other, the measurement of slowly growing instabilities in
numerical simulations also has a number of technical difficulties.
In particular, measuring the growth of these very weak instabilities
above the background noise requires a sufficient integration time,
and in *N*-body simulations as well as in actual physical systems,
the time-scale required for the growth of such instabilities can
overlap with the secular time-scale of two-body relaxation in the
cluster. As a consequence, one might miss the growth of an insta-
bility, because the system has already secularly evolved to a stable
configuration.

Yet, even if we only focus on reasonably large growth rates (say
η > 0.05), the results of Fig. 10 show that the stability of spherical
clusters has at least one additional degree of freedom compared to
the single one used by Ostriker & Peebles (1973). As discussed in
Section 5.1.2, rotation in the innermost part of the cluster seems to
have a critical influence on the increase in growth rate, compared to
the non-rotating case. Fig. 11 shows that a combination of *T/W* and the
mean azimuthal velocity at *R* = 0.2 makes a better stability criterion
than the original *T/W*, even if the separation between stable and
unstable systems (or, rather, slowly and rapidly growing instabilities)
in that space is not yet perfect.

Let us remark in closing that the 2D space of global rotation
(through *T/W*) and azimuthal velocity at a given radius can be mapped
into the space of global rotation and shear: a cluster with high
|**v**<sub>φ</sub>(0.2)| displays low shear if *T/W* is also high, and high shear if
*T/W* is low; a cluster with low |**v**<sub>φ</sub>(0.2)| displays high shear if *T/W* is
high, and low shear if *T/W* is low.

5.2.2 Non-rotating models with nearly circular motions

We take the view that the question of the stability of a non-rotating
equilibrium with a high degree of tangential anisotropy is still
unsettled. On the one hand, Polyachenko (1987) uses the WKB
approximation to show that such spheres are linearly stable. However,
this result only probes the case of quasi-local disturbances, while we
measure large-scale patterns (as illustrated by Fig. 13), which may be
beyond the scope of the WKB approximation. On the other hand, the
possibility of large-scale instabilities in such spheres is supported
by Palmer et al. (1989), based on complementary analytical and
numerical results. A consequence of this finding, if confirmed, would
be that sufficiently tangentially anisotropic clusters seem to display
unstable behaviour, *whatever their degree of rotation or shear*. But
we caution that this intriguing regime deserves a careful analysis that
goes beyond the scope of the current paper, and we hope to present
a contribution dedicated to this topic in the future.

5.3 Conclusions

This paper focuses on the collisionless stability of rotating, spherical stellar systems, using a mix of *N*-body simulations, the matrix method for the evolution of perturbations, and some physical ideas. Our starting point is a set of anisotropic Plummer models (Dejonghe 1987), whose anisotropy is parametrized by a parameter *q*; our focus is on models with *q* < 0, whose velocity distribution is tangentially biased. 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 integrals of motion, and parametrized in various ways by a parameter *γ*. Instances of these models were generated with the publicly available software PlummerPlus, and their evolution studied using the public version of gyrfalcON.<fnref n="9" /> Supplementary studies used NBODY6 (Nitadori & Aarseth 2012). The matrix method, described for example in Rozier et al. (2019), was also used to search for unstable modes.

Our conclusions are as follows:

(i) For values of *q* < 0 and 0 < *α* ≤ 1, the models are rather generally unstable to the formation of a rotating ‘bar’ perpendicular to the axis of rotation.

(ii) In models in which the anisotropy parameter $q = -6$, and a fraction $\alpha$ of retrograde stars have the rotation reversed, the growth rate of the mode is approximately proportional to $\alpha^6$ (Fig. 5), at least when $\alpha$ is large enough that the growth rate can be measured reliably. In extremely tangentially anisotropic models, the pattern speed diminishes as $\alpha$ 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 $\gamma_1$ of the minimum energy have their rotation reversed, the growth rate reaches roughly its maximum value if $\gamma_1 \gtrsim 0.5$, i.e. if the Lynden-Bell demon is applied to at least the most bound half of the stars (Fig. 6).

(iv) In models in which stars below a given energy have the Lynden-Bell demon applied in a prograde sense and for the remainder it is applied in a retrograde sense, the minimum growth rate occurs when the proportions of retrograde and prograde stars in the most bound 40 per cent of the system are nearly the same (i.e. the case $\gamma_2 = 0.2$ in Fig. 9).

(v) If the Lynden-Bell demon is applied to only the stars with
orbital inclination above (or below) some limit, the growth rate

<footnote n="9">⁹https://teuben.github.io/nemo/man.html/gyrfalcON.1.html</footnote>
continues to increase as the limit increases to include polar orbits (or decreases to include equatorial orbits, respectively) (Fig. 12).

(vi) The bar produced by the instability is detectable only within the cylinder $R \lesssim 1.5$, $|z| < 0.9$ in Hénon units (Section 3.4.1).

(vii) Though the dynamics is affected by several other resonances beyond the familiar Lindblad and corotation resonances, no individual resonance appears to play an essential role (Table 2).

(viii) The pattern speed of the bar, which varies surprisingly little from one model to another, is such as to avoid the ILR everywhere, and the corotation resonance in the cylindrical region where the bar is found (Figs 14 and 15).

(ix) The ILR is the single most important resonance in the bar instability, even though there are no particles at the ILR (Fig. 16 and Table 2). Other significant resonances are the OLR, the corotation resonance, and two resonances that we refer to as 'tumbling' and 'epicyclic' resonances (Table 2 and Figs 17 and 18).

(x) It cannot be shown, from the results in this paper, that any of the surveyed models is stable against bar formation. If, instead of searching for a stability criterion, one searches for a criterion that determines which models may exhibit rapid instability, we show that, in addition to the global rotational kinetic energy, it is helpful to include a measure of the rotation speed at a radius inside the virial radius of the system (Fig. 11).

## ACKNOWLEDGEMENTS

We thank an anonymous referee for a number of helpful comments, including pertinent references. SR thanks Christophe Pichon and Jean-Baptiste Fouvry for many interesting discussions, and all authors are grateful to them for commenting on an early draft. This work has made use of the Horizon Cluster hosted by Institut d'Astrophysique de Paris. All authors thank Stephane Rouberol for running smoothly this cluster. We are grateful to S. Law and his IT staff in the Mathematics Department at Edinburgh University, for very helpful management of the GPU hardware. Many simulations were also carried out with the ECDF (Edinburgh Compute and Data Facility) cluster at the University. All authors thank the Leverhulme Trust for financial support from grant RPG-2015-408. ALV acknowledges support from a UKRI (United Kingdom Research and Innovation) Future Leaders Fellowship (MR/S018859/1) and is grateful to the Department of Astronomy at the University of Tokyo, where part of this work was conducted, thanks to a JSPS (Japan Society for the Promotion of Science) International Fellowship with Grant-in-Aid (KAKENHI-18F18787).

## DATA AVAILABILITY

Data underlying this article will be shared on reasonable request to the authors via d.c.heggie@ed.ac.uk. Data related to the initial conditions may be reproduced via the publicly available software PlummerPlus.

## REFERENCES

Alimi J.-M., Perez J., Serna A., 1999, MNRAS, 305, 859

Allen A. J., Palmer P. L., Papaloizou J. C. B., 1992, MNRAS, 256, 695

Anderson J., King I. R., 2003, AJ, 126, 772

Antonov V. A., 1987, in de Zeeuw P. T., ed., Proc. IAU Symp. 127, Structure and Dynamics of Elliptical Galaxies. Kluwer, Dordrecht, p. 531

Aoki S., Noguchi M., Iye M., 1979, PASJ, 31, 737

Athanassoula E., 2008, MNRAS, 390, L69

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

Bellini A. et al., 2017, ApJ, 844, 167

Berman R. H., Mark J. W.-K., 1979, A&A, 77, 31

Bertin G., Lin C. C., 1996, Spiral Structure in Galaxies: A Density Wave Theory. MIT Press, Cambridge, MA

Bertin G., Pegoraro F., Rubini F., Vesperini E., 1994, ApJ, 434, 94

Bianchini P. et al., 2018, MNRAS, 481, 2125

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

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

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

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

Christodoulou D. M., Shlosman I., Tohline J. E., 1995, ApJ, 443, 551

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

Dejonghe H., 1987, MNRAS, 224, 13

Dejonghe H., Merritt D., 1988, ApJ, 328, 93

Doremus J.-P., Feix M. R., Baumann G., 1971, Phys. Rev. Lett., 26, 725

Efstathiou G., Lake G., Negroponte J., 1982, MNRAS, 199, 1069

Evans N. W., Read J. C. A., 1998, MNRAS, 300, 106

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

Fujiwara T., 1983, PASJ, 35, 547

Gaia Collaboration, 2018, A&A, 616, A12

Gillon D., 1980, CRASB, 290, 545

Hachisu I., 1986, ApJS, 61, 479

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

Heggie D. C., Breen P. G., Varri A. L., 2020, MNRAS, 492, 6019

Hénon M., 1973, A&A, 24, 229

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

Hohl F., 1971, ApJ, 168, 343

Hohl F., 1976, AJ, 81, 30

Jalali M. A., Hunter C., 2005, ApJ, 630, 804

James R. A., 1964, ApJ, 140, 552

Jindal A., Webb J. J., Bovy J., 2019, MNRAS, 487, 3693

Kalnajs A. J., 1972, ApJ, 175, 63

Kalnajs A. J., 1977, ApJ, 212, 637

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

Karino S., Eriguchi Y., 2003, ApJ, 592, 1119

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

Lai D., Rasio F. A., Shapiro S. L., 1993, ApJS, 88, 205

Lin C. C., 1955, The Theory of Hydrodynamic Stability. Cambridge Univ. Press, Cambridge

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

Lynden-Bell D., 1962, MNRAS, 123, 447

Lynden-Bell D., 1969, MNRAS, 144, 189

Lynden-Bell D., Kalnajs A. J., 1972, MNRAS, 157, 1

Lynden-Bell D., Sanitt N., 1969, MNRAS, 143, 167

Mark J. W.-K., 1977, ApJ, 212, 645

Merritt D., 1999, PASP, 111, 129

Merritt D., Aguilar L. A., 1985, MNRAS, 217, 787

Meza A., 2002, A&A, 395, 25

Mikhálovskií A. B., Fridman A. M., Epel'Baum Y. G., 1971, Sov. J. Exp. Theor. Phys., 32, 878

Milone A. P., Marino A. F., Mastrobuono-Battisti A., Lagioia E. P., 2018, MNRAS, 479, 5005

Morozov A. G., Polyachenko V. L., Shukhman I. G., 1974, SvA, 18, 44

New K. C., Centrella J. M., Tohline J. E., 2000, Phys. Rev. D, 62, 064019

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

Orr W. McF., 1907, Proc. Roy. Irish Acad. A, 27, 9

Ostriker J. P., Mark J. W.-K., 1968, ApJ, 151, 1075

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

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

Palchik M. I., Patashinskii A. Z., Pinus V. K., 1974, PMTF, 63

Palmer P. L., 1994, Stability of Collisionless Stellar Systems. Kluwer, Dordrecht

Palmer P. L., Papaloizou J., 1987, MNRAS, 224, 1043

Palmer P. L., Papaloizou J., Allen A. J., 1989, MNRAS, 238, 1281

Palmer P. L., Papaloizou J., Allen A. J., 1990, MNRAS, 243, 282

Perez J., Aly J.-J., 1996, MNRAS, 280, 689
Perez J., Alimi J.-M., Aly J.-J., Scholl H., 1996, MNRAS, 280, 700

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

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

Polyachenko V. L., 1987, in de Zeeuw P. T., ed, Proc. IAU Symp. 127, Structure and Dynamics of Elliptical Galaxies. Kluwer, Dordrech, p. 301

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

Polyachenko E. V., 2005, MNRAS, 357, 559

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

Polyachenko E. V., Polyachenko V. L., Shukhman I. G., 2007, MNRAS, 379, 573

Rozier S., Fouvry J.-B., Breen P. G., Varri A. L., Pichon C., Heggie D. C., 2019, MNRAS, 487, 711

Saha P., 1991, MNRAS, 248, 494

Saijo M., Yoshida S., 2006, MNRAS, 368, 1429

Saijo M., Yoshida S., 2016, Phys. Rev. D, 94, 084032

Saijo M., Baumgarte T. W., Shapiro S. L., 2003, ApJ, 595, 352

Sellwood J. A., 1981, A&A, 99, 362

Sellwood J. A., 1997, in Clarke D. A., West M. J., eds, ASP Conf. Ser. Vol. 123, Computational Astrophysics. Astron. Soc. Pac., San Francisco, p. 215

Sellwood J. A., Valluri M., 1997, MNRAS, 287, 124

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

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

Stoeckly R., 1965, ApJ, 142, 208

Synakh V. S., Fridman A. M., Shukhman I. G., 1971, Dokl. Akad. Nauk SSSR, 201, 827

Tassoul J.-L., 1978, Theory of Rotating Stars. Princeton Univ. Press, Princeton, NJ

Tohline J. E., Durisen R. H., McCollough M., 1985, ApJ, 298, 220

Toomre A., 1981, in Fall S. M., Lynden-Bell D., eds, The Structure and Evolution of Normal Galaxies. Cambridge Univ. Press, Cambridge, p. 111

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

Vandervoort P. O., 1982, ApJ, 256, L41

Vandervoort P. O., 1991, ApJ, 377, 49

Vandervoort P. O., Welty D. E., 1982, ApJ, 263, 654

Watts A. L., Andersson N., Beyer H., Schutz B. F., 2003, MNRAS, 342, 1156

Watts A. L., Andersson N., Williams R. L., 2004, MNRAS, 350, 927

Watts A. L., Andersson N., Jones D. I., 2005, ApJ, 618, L37

Weinberg M. D., 1994, ApJ, 421, 481

Yoshida S., Saijo M., 2017, MNRAS, 466, 600

Zuang T. A., Hohl F., 1978, ApJ, 226, 521

prevented the determination of the pattern speed by the same method. For the pattern speed, data from several runs could not be simply combined, and an alternative method was used, viz. calculation of the autocorrelation of the time series of $\arg(C_2)$. Averaging results over several runs (as mentioned above) gave a reliable autocorrelation with an estimate of the standard error. The form of the average autocorrelation function is a clear decaying oscillation. The location of its first maximum was determined by applying a quadratic fit to the section on either side of the maximum, as far as the neighbouring minima. This method does not allow determination of the sense of pattern rotation. This is therefore undetermined for $\alpha = 0.2$, but is presumed to be prograde. For $\alpha = 0$, it may be either, or possibly a mixture (actually surface-harmonic expansion shows that there are five modes with $l = 2$, but three are filtered out by the definition of $C_2$).

## APPENDIX B: RESONANCES IN INCLINED EPICYCLIC MOTIONS PERTURBED BY A ROTATING BAR

In the context of disc dynamics, Binney & Tremaine (2008, section 3.3.3) use the expression $\Phi_1(R, \varphi) = \Phi_b(R)\cos[2(\varphi - \Omega_b t)]$ as a model for the potential of a planar bar, where $R, \varphi$ are plane polar coordinates in the orbital plane, $\Omega_b$ is the pattern speed, and $\Phi_b$ specifies the radial dependence of the potential. The obvious extension to three dimensions is $\Phi_1(r, \theta, \varphi) = \Phi_b(r)\sin^2\theta\cos[2(\varphi - \Omega_b t)]$ in spherical coordinates, since the angle-dependence is that of an $l = 2$ spherical harmonic.

We consider the first-order perturbing potential of a particle which, at lowest order, is in epicyclic motion on a plane inclined at inclination $i$ to the x, y plane ($\theta = \pi/2$), and passes its ascending node at longitude $\Omega$ at time $t = 0$.<fnref n="10" /> In cylindrical polar coordinates $(R, \varphi, z)$ based on the orbital plane, with the angular coordinate $\varphi = 0$ at the node, it is easy (but a bit tedious) to show that

$$ \Phi_1(R, \varphi, z=0) = \frac{1}{4} \Phi_b(R) \left\{ 2 \sin^2 i \cos[2(\Omega - \Omega_b t)] + (1+\cos i)^2 \cos([2(\Omega - \Omega_b t + \varphi)] + (1-\cos i)^2 \cos([2(\Omega - \Omega_b t - \varphi)]) \right\}. \quad (B1) $$

Now we make use of the epicyclic approximation, by which (Binney & Tremaine 2008, equations 3.91 and 3.93a)

$$ R = R_0 + X \cos \theta_1 \quad (B2) $$

$$ \varphi = \theta_2 - \gamma X \sin \theta_1, \quad (B3) $$

where $(R_0, \theta_2)$ are the cylindrical radius and longitude of the guiding centre, $X$ is the radial amplitude of the epicyclic motion, $\theta_1$ is the phase of the epicycle, and $\gamma$ is a certain function of $R_0$. Substituting equations (B2) and (B3) into equation (B1) and expanding to lowest order in $X$, we can readily locate all trigonometrical terms of the form

$$ \left\{ \begin{matrix} \cos \\ \sin \end{matrix} \right\} (2(\Omega - \Omega_b t) + n_1 \theta_1 + n_2 \theta_2), \quad (B4) $$

where $n_1, n_2$ are (small) integers. This includes the pair $(0,0)$, which arises from the first term in equation (B1). Using the fact that $\dot{\theta}_{1,2} = \Omega_{1,2}$ (the radial and azimuthal frequencies), these terms yield the resonances listed in Table 1, including the stated $i$-dependence and the values in the last column.

## APPENDIX A: DETERMINATION OF TIME-SCALES IN THE MODELS OF SECTION 3.4.2

This appendix describes details of the way in which data for Fig. 14 were derived from results of N-body simulations.

The analysis of $C_2$ (equation 2) began with the growth rate. For large $\alpha$, $|C_m|$ was found to saturate at some time $t_{max}$. For these models ($\alpha \ge 0.6$), the growth rate was obtained by fitting a logistic curve (which is nearly exponential at early times) in this time interval. For smaller $\alpha$, where the growth of a bar (if present) did not saturate by time $t = 200$, a simple exponential fit was used. In these cases, however, data from several models were combined: eight for $\alpha = 0.4$ and 10 each for $\alpha = 0, 0.2$. In all cases, error bars represent the $1\sigma$ asymptotic confidence interval.

The pattern speed was estimated by examining the time series of $\arg(C_2)$. For a rotating bar, this exhibits a sawtooth wave form, the angular speed of each segment being twice the pattern speed. For large $\alpha$ it was evaluated by simply counting the number of periods within the selected time interval (see above), and the motion could be easily seen to be prograde. For $\alpha = 0.4$, the expected prograde pattern was still visible in places, but with irregularities which

<footnote n="10">¹⁰In the rest of the paper, $\Omega$ means the circular angular velocity.</footnote>
APPENDIX C: APPLICATION OF THE
RESTRICTED MATRIX TO THE RADIAL ORBIT
INSTABILITY

As described in Binney & Tremaine (2008), spherical systems biased towards radial orbits are subject to the so-called radial orbit instability. This instability emerges through a Jeans process, by the clustering of radially elongated orbits around a bar-like overdensity. More precisely, the torque applied by a weak bar on a radially elongated orbit acts to modify the apsidal precession rate of the orbit (i.e. $\Omega_2 - \Omega_1/2$), in such a way that the angle between the bar and the ellipse tends to decrease. This process results in a bar-like instability, with a pattern speed in the range of the radial orbits' precession rates. In this process, the orbits around the ILR are particularly critical, since they form the main body of the instability. In a rotating system, Palmer et al. (1989) show that the tumbling process also participates in the destabilization of the radial orbit instability. Here, we apply the restricted matrix method to the case of a radially anisotropic, rotating Plummer sphere, which gives a well-constrained test-case for its application.

Table C1 reports the values of the growth rate and oscillation frequency found by the restricted matrix method in the $q = 2, \alpha = 1$ system, when a variety of combinations of resonance vectors are tested. The computation of reference includes all relevant resonance vectors with $|n_1|, |n_2| \le 2$, corresponding to $\mathcal{I}_0 = \{(0, 0), (1, 0), (2, 0), (-2, 2), (-1, 2), (0, 2), (1, 2), (2, 2)\}$. It appears that turning off terms corresponding either to the ILR

<table><caption>**Table C1.** Values of the growth rate and pattern speed found through the truncated response matrix method applied to the $q = 2, \alpha = 1$ radially biased system.</caption><thead><tr><th><i>I</i></th><th><i>η</i></th><th><i>Ω</i><sub>p</sub></th></tr></thead><tbody><tr><td>Reference</td><td>0.088</td><td>0.073</td></tr><tr><td>(-1, 2), (0, 0)</td><td>0.088</td><td>0.066</td></tr><tr><td><i>I</i><sub>0</sub> \ (-1, 2)</td><td>&lt;2 × 10<sup>-4</sup></td><td>—</td></tr><tr><td><i>I</i><sub>0</sub> \ (0, 0)</td><td>0.033</td><td>0.027</td></tr><tr><td>(-1, 2)</td><td>0.018</td><td>0.024</td></tr><tr><td>(0, 0)</td><td>&lt;2 × 10<sup>-4</sup></td><td>—</td></tr></tbody></table>

or to the tumbling process from the full matrix lowers the growth rate and modifies the oscillation frequency of the instability, while keeping only those two terms together preserves both its growth rate and (nearly) its oscillation frequency. These results show that the destabilization of radially anisotropic rotating systems is guaranteed through the interaction of the radial orbit and tumbling processes only. In that case, the identification of these processes is made easy by the largely dominant role played by the corresponding terms in the response matrix: the second line of Table C1 shows that all other resonance vectors can be neglected w.r.t. (-1, 2), (0, 0). Section 4.2 shows that tangentially anisotropic systems are more complex in that respect.

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